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SUMMARY 

The  principal  aim  of  the  thesis  is  to  develop,  and  shew  how  to  solve, 
equations  describing  the  flexural  notions  caused  in  a  ship  by  the  pulsa¬ 
tions  of  the  bubble  of  gaseous  products  from  a  nearby  underwater  explosion. 
Equations  for  this  purpose  were  developed  in  the  USA  15-20  years  ago  and 
the  thesis  presents  an  alternative,  finite  element  formulation  which  is 
more  general  and  which  leads  to  equations  in  matrix  form.  In  this  form 
the  equations  are  much  more  suitable  for  evaluation  hy  computer.  The 
method  of  solution  of  these  equations  of  elastic  motion,  in  terms  of 
normal  mode  shapes,  is  described. 

The  approximations  involved  in  the  hydrodynamics  of  the  bubble-ship  inter¬ 
action  are  examined  in  some  detail  and  the  limits  of  application  of  the 
theory  are  illustrated.  The  principal  approximation  in  the  hydrodynamics 
is  the  use  of  strip  theory  since  this  cannot  allow  for  either  the  longitud¬ 
inal  flow  which  occurs  when  the  bubble  is  close  to  the  ship,  or  f or  the 
flow  around  the  end  of  the  ship  when  the  bubble  is  in  such  a  region.  The 
accuracy  of  strip  theory  is  examined  in  both  these  cases.  Ar.  alternative 
to  strip  theory  for  describing  the  flow  around  a  vibrating  ship  is  presented. 
This  approximates  the  full  three-dimensional  flow  field  around  the  ship  by 
means  of  a  set  of  line  distributions  of  vertically  oriented  dipoles  at  the 
waterplane-ship  centreline  axis.  The  technique  could  also  be  of  value  for 
the  determination  of  added  water  mass  around  a  vibrating  ship  in  more  usual 
vibration  problems. 

The  original  American  derivation  of  the  equations  assumed  that  the  pulsating 

bubble  was  non-migrating.  Suitable  equations  are  developed  to  describe  the 

/ 

upward  migration  of  the  bubble  due  to  gravity  and  the  effect  of  this  migra¬ 
tion  on  the  motion  of  the  ship. 

For  severe  flexural  motions  of  the  ship  some  plastic  yielding  or  even  buck¬ 
ling  could  occur  and  a  direct  solution  of  the  original  finite  element 
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equations,  modified  to  include  such  effects,  is  presented.  This  shows  thi 
way  in  which  they  can  change  the  motion  of  the  ship. 

The  entire  thesis  is  a  theoretical  exercise  although  the  results  for  the 
non-migrating  cubble  case  are  compared  with  a  few  model  test  results 
published  in  connection  with  the  original  American  work. 
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Introduction 


When  a  fairly  large  explosive  charge  is  detonated  underwater 
close  to  a  ship  the  shock  wave  it  transmits  to  the  water  can  do 
considerable  damage  to  the  ship.  If  close  enough  the  hull  plating  say 
rupture  and  for  charges  farther  away  there  may  still  be  considerable 
deformation  of  the  structure  nearest  the  charge.  The  shock  wave  also 
is  partially  transmitted  into  the  structure  and  can  cause  farther 
shock  damage  to  internally  mounted  equipment.  These  effects  follow 
close  behind  the  detonation,  typically  witnin  times  of  the  order  of 
15  to  20  milliseconds,  and  are  usually  confined  to  the  region  of  the 
ship  closest  to  the  explosion. 

The  explosion  may  however  also  excite  a  very  different  kind  cf 
response.  Ships  near  underwater  explosions  often  sheke  and  bend 
overall  in  violent  flexural  vibrations.  Such  vibrations  are  usually 
confined  largely  to  a  vertical  plane  and,  although  very  violent,  are 
of  much  lower  frequency  (typically  1-5  Hz)  than  the  shock  motions 
accompanying  the  local  damage.  The  amplitudes  can  be  veiy  large  and 
can  cause  considerable  structural  damage  to  the  ship.  Unlike  the 
direct  shock  damage,  this  damage  can  occur  at  parts  of  the  ship 
quite  remote  from  the  site  of  the  explosion.  These  bending  motions 
are  commonly  termed  'whipping'  notions.  They  also  occur  in  severe 
storms  if  the  bow  of  a  ship  emerges  from  the  water  and  then  slams 
back  onto  it.  They  constitute  a  superposition  of  the  well  known  low 
frequency  vertical  bending  vibration  modes  of  the  ship,  excited  to 
unusually  large  amplitudes.  It  is  the  aim  of  this  thesis  to  investi¬ 
gate  the  nature  of  the  interaction  between  the  explosion  and  the  ship 
which  causes  these  whipping  motions,  and  to  develop  suitable 
equations  from  which  whipping  motions  can  be  calculated. 
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It  is  already  known  from  earlier  work  principally  by  G  Chertock 
(and  described  briefly  in  Section  B)  that  such  notions  are 
connected  closely  to  the  notion  of  the  water  around  the  bubble  of 
gaseous  combustion  products  fron  the  explosion.  This  bubble  i3  well 
known  to  pulsate  vigorously  and  the  period  of  the  pulsation  is  often 
very  close  to  the  low  frequency  vibration  nodes  of  ships.  Chertock*s 
work  showed  that  the  whipping  is  primarily  due  to  the  coupling  between 
the  vibration  nodes  and  the  pulsating  flow.  This  thesis  extends  the 
earlier  work,  investigates  the  hydrodynamic  interaction  between  the 
explosion  and  the  ship  in  wore  detail  and,  by  using  a  different 
initial  approach,  provides  a  set  of  equations  for  calculating  the 
ship  response  which  show  more  clearly  the  nature  of  the  forces 
acting  on  the  ship  and  are  also  nore  readily  solved  by  standard 
matrix  methods  on  a  computer. 

The  mechanics  of  the  interaction  involves  several  standard 
fields  of  study;  the  elastic  nature  of  the  ship,  the  tydrodynanics 
of  vibration  and  the  hydrodynamics  of  underwater  explosions.  The 
thesis  attempts  to  tie  together  only  the  simplest  ideas  from  each 
and  to  extend  standard  methods  in  each  only  where  they  appeared  to 
be  deficient  from  a  whipping  viewpoint:  for  example  the  strip  theoiy 
of  ship  vibration  and  the  equations  for  a  pulsating  migrating 
bubble.  In  consequence  the  bibliographies  given  in  each  field  are 
rather  basic;  they  are  aimed  simply  at  providing  sufficient  back~ 
ground  information  to  apply  adequately  the  methods  developed.  The 
first  part  0f  this  thesis,  sections  A  to  K,  is  largely  based  on 
existing  methods,  although  the  methods  have  not  previously  been 
used  for  this  particular  application, and  is  concerned  largc’y  with 
establishing  the  equivalence  of  the  present  methods  to  the  earlier  work. 
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Section  D  extends  strip  theory  to  allow  for  the  explosion  induced  motion 
of  the  surrounding  water.  Section  G  includes  a  simple  approach  to  making 
an  allowance  for  the  effect  of  smal  1  bilge  keels  on  the  added  water  mass. 
Section  I  describes  a  large  computer  program  written  to  evaluate  thi 
equations  developed  in  sections  C  to  H,  Most  of  the  new  material  in  the 
thesis  is  contained  in  sections  L  to  0  and  in  11.  Section  L  describes  an 
alternative  to  the  customary  ship  theory  for  the  flow  around  the  ship.  « 
Section  M  investigates  the  limits  of  application  of  the  approximate  analysis 
by  comparing  the  forces  predicted  by  the  approximate  method  with  exact 
solutions  for  very  simple  cases.  Section  N  develops  an  approximate  me  thou 
for  making  an  allowance  for  the  effect  of  shock  wave  damage  to  the  target. 
Section  0  investigates  the  effect  of  gravity  on  the  explosion  bubble  and 
on  the  forces  the  bubble  produces  on  the  ship.  Section  R  describes  how  the 
basic  equations  for  the  ship  motion  can  be  modified  to  allow  for  plastic 
yielding  and  buckling  of  the  ship.  Two  further  major  programs  were 
developed  to  evaluate  the  ship  response  (including  allowances  for  gravity 
migration  of  the  bubble  and  plastic  deformation  of  the  target)  and  are 
described  in  sections  P  and  R.  Several  minor  computer  programs  were 
developed  to  evaluate  the  exact  solutions  for  comparison  with  the  approxi¬ 
mate  analysis. 

To  make  each  section  of  the  thesis  as  self-contained  as  possible  the  refer¬ 
ences  and  figures  for  each  section  are  grouped  at  the  end  of  the  section. 


8  Background 


Although  a  coiwiderable  Mount  of  work  has  been  carried  out  on 
the  vibration  and  «lacaing  aspects  of  ship  whipping*  very  little  work 
has  been  published  on  the  explosion  induced  whipping  notions.  What 
work  has  been  published  is  due  entirely  to  a  single  author*  G  Chertock. 
His  work  is  described  in  references  fitt]  to  [B3l but,  as  it  appears  to  be 
the  sole  predecessor  of  the  work  here  it  is  worth  describing  it  in 
s<mb  detail. 

The  ship  target  is  assuaed  to  be  equivalent  to  an  elastic  free- 
free  bean  whose  transverse  displaceaents  can  therefore  be  represented 

by 

y(x,t)  =  |qi(t)¥.(x)  (Bi) 

where  x  is  the  longitudinal  distance  along  the  axis  of  the  ship* 
y  the  vertical  displacement  of  the  cross-section  at  x  (assuaed  rigid)* 
¥^(x)  the  ith  normal  mode  shape  for  transverse  vibration  in  vacwo 
and  qj  is  a  generalised  co-ordinate  for  the  i*'*1  node  (e.g.  the 
amplitude  at  the  bow  or, alternatively* the  root  mean  square  displace¬ 
ment).  The  normal  modes  for  such  motions  are  orthogonal  with  the 
mass  per  unit  length,  m(x),  as  a  weighting  function. 

/^m(x)¥.(x)¥j(x)dx  =  H.6^  (B2) 

where  1L  is  the  generalised  mass  associated  with  the  i^  mode. 

By  using  the  orthogonality  condition  the  equations  of  vertical 
transverse  motion  due  to  an  arbitrary  distribution  of  pressure  over 
the  surface  of  the  ship  can  be  written 

*Vqi  +  <xiqi^  =  Qi  =  “  /lsP(s,t)¥.cos3  da  i  =  0,1,2* w  (B3) 
where  is  the  circular  frequency  for  mode  i,  the  generalised 
force  and  da  an  element  of  area  of  tne  surface  of  the  ship.  The 
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angle  p  is  the  angle  between  the  vertical  and  the  normal  to  the 

surface  and  s  refers  to  a  point  on  the  ship  surface. 

By  neglecting  the  Bernoulli  pressures,  -'^pu  ,  compared  with 

tiie  unsteady  forces  p|^  in  the  usual  Bernoulli  pressure  equation 

for  incompressible  flow,  the  fluid  pressure  on  the  ship  can  be  expressed 

as  a  linear  sum  of  the  pressure  due  to  the  explosion  induced  flu'.d  flow 

on  the  structure,  assumed  rigid,  and  the  pressures  due  to  the  motion 

of  the  target  in  its  various  modes,  so  that 

p(s,t)  =  P0Q(s,t)  ♦  Sp.(s,t)  (B4) 

where  p  is  due  to  the  explosion  and  p.  to  the  ith  mode. 

With  this  assumption  of  linearity,  induced  modal  mass  coefficient 

can  be  defined  by 

L.  .  =  -  [//  p.(s,tn;  cospdoj/q  (t)  (B5) 

xj  5  J  X  A 

and  will  be  independent  of  time.  Substituted  into  equation  (B3) 
this  gives 

(K.  +  L^)  (q  i  +  co^qj)  =  -  /Jp^^cos^dcr  -  U,.q.;  i  =  0,1,2,... 


1 .  +  L .  . 
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In  these  equations  the  elastic  modes  are  clearly  coupled  to  each 
other  by  the  entrained  water  around  the  ship  through  the  coefficients 
i  £  j.  By  a  suitable  transformation  of  the  normal  mode  shapes 
^  and  co-ordinates  it  was  shown  Uiat  the  equations  can  be 
uncoupled  into 

•  M*. ( q /  +  w  f 2q .')  =  -  / /p  T.'cos^do  (B7) 

i  i  i  V  oo  l  r 

The  procedure  however  is  extremely  difficult  to  carry  out  in  practice 
due  to  the  necessity  to  evaluate  the  cross-coupling  coefficients. 
These  can  only  be  determined  exactly  by  the  solutions  of  .a  series 
of  integral  equations  where  the  integrals  are  over  the  whole  surface 
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of  the  ship  hull.  In  practice,  for  physical  reasons  the  L.  ,'s 

^  J 

are  expected  to  be  fairly  snail  if  the  ship  is  long  and  slender  or 
if  there  is  sufficient  fore-and-aft  symmetry.  Chertock  avoided 
determination  of  the  I„.  .  by  assuming  that  ships  are  'proportional 
bodies' 

i.e.  by  his  definition,  cos(5  ds  =  a^pmtx'^Cx) 
where  the  integral  is  a  line  integral  taken  round  a  ship  cross- 
section,  p  is  the  density  of  water,  and  the  relation  applies  at  all 
ship  cross-sections,  i.e.  for  all  i,a^  is  assumed,  for  proportional 
bodies,  to  be  a  constant,  independent  of  x.  With  this  definition, 
by  multiplying  by  ¥ .  and  integrating  over  the  length  of  the  ship. 


L.  .  = -^-//p -(3,t)¥.cospdcr  =  -  .  <jx£p .  cospds  = - — J^mT.V.dx 

1J  ••  J  *  •»  t)  1  J  •  •  O  J  X 


3  a  ,P 

i.e.  L.  .  =  -  — M.6.  . 
U  *  i  iJ 

j 


"j 


qj 


so  that  the  cross-coupling  terms  all  vanish  and  the  mode  shapes  in 
air  and  water  are  identical,  although  the  frequencies  will  change 
due  to  the  non-zero  entrained  water  masses 

IX 

For  the  term  p^o,  representing  the  pressure  on  the  ship,  if 

it  were  rigiu,  due  to  the  flow  around  the  explosion  bubble,  the  well 

known  result uB/,]tliat  the  flow  around  a  non-migrating  bubble  is 

equivalent  to  that  of  a  simple  source  was  used.  Since  the  rate 

♦ 

of  change  of  the  bubble  volume  is  V(t),  the  strength  of  the  source 
is  V/4n  and  its  potential  is 


.  -  _L 

'o  4 3tr 

The  pressure  in  the  absence  of  the  target  is  therefore 
•  • 

po  ~  P^o  “  4^r 

ami  the  integral  -//poo¥  jCosJido  will  be  representable  in 

•  • 

pV 


(B3) 

the  form 
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whore  p  is  a  factor  depending  only  on  the  charge  position  and  the 
ship  and  mode  shapes.  For  a  cliarge  sufficiently  far  from  the  ship 
the  factor  p  was  shown  to  be  given  approximately  by 


-  ,  uL. .  7. 

■i s  (pa  *  ttH"1* 


(B9) 
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where  £  is  tuc  lieight  of  the  axis  of  the  ship  above  the  charge  centre 
and  r  is  the  standoff  of  the  particular  point  on  the  ship  axis  from 
the  charge  centre.  The  final  modal  equation  of  motion  is  therefore 


(Mi +  Lii)(iii  ^iV 


(B10) 


To  determine  the  response  of  the  ship  then,  equation  (BIO) 

is  to  be  solved  for  each  mode,  using  frec-field  values  for  the 

•  • 

explosion  bubble  volume  acceleration,  V,  and  the  modal  components 
may  then  be  summed  to  give  the  displacement  function  y(x,t).  The 
principal  assumptions  made  in  the  derivation  are: 

a.  The  charff  is  far  enough  from  the  ship  for  the 
explosion  bubble  motion  not  to  be  affected 
by  the  ship. 

b.  The  ship  motions  remain  completely  elastic. 

c.  The  ship  is  a  'proportional  body',  i.e.  there 

is  no  fluid  coupling  between  the  'in  vacuo'  elastic 
vibration  modes  when  the  ship  is  in  the  water. 

d.  The  charge  is  at  a  distance  from  the  ship  large  com¬ 
pared  with  the  ship  cross-section,  so  that  the 
approximation  (B9)  applies. 

e.  There  is  no  significant  migration  of  the  bubble 
during  the  motion. 

f.  Only  incompressible  water  motions  are  of  importance. 

g.  Structural  and  hydrodynamic  damping  are  negligible. 

To  check  the  accuracy  of  his  equations,  Chertock  carried  out 
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two  scries  of  experiments  on  model  targets.  The  first  series 
used  a  1.2  gn  charge  under  a  floating  rectangular  box  10  ft  long 
by  9  ins  wide  by  6  ins  deep.  This  target  was  slender,  symmetric 
and  nearly  uniform  and  the  mode  shapes  in  air  and  water  were 
indistinguishable  so  that  the  target  certainly  satisfied  the 
requirement  v,c).  The  charge  depth  of  13  ins  was  carefully  selected 
so  that  the  normal  repulsion  from  a  free  surface  exactly  balanced 
the  upward  buoyancy  migration  due  to  gravity  so  that  requirement  (e) 
was  also  satisfied.  At  the  given  depth  the  maximum  bubble  radius 
is  about  6  ins  so  tliat  (a)  is  satisfied  and  the  geometry  is  clearly 
such  that  (d)  is  not  really  violated  although  it  might  be  thought 
near  the  limits  of  acceptability  since  the  cross-sectional  dimensions 
are  nearly  comparable  to  the  standoff.  For  the  fundamental,  2-node 
whipping  mode  the  predicted  results  agreed  very  well  with  the  experi¬ 
mental  (Reference  B3)  but  for  the  J-node  mode  the  measured  amplitudes  were 
about  30  less  than  those  predicted.  For  the  2-node  mode  the  varia¬ 
tion  of  whi  ping  amplitude  with  the  longitudinal  location  of  the 
charge  agreed  very  weli  with  experiment. 

The  second  series  of  tests  used  a  submerged  cylindrical  model 
8  ft  long  by  9  ins  in  diameter  and  two  charge  sizes,  I .2  gms  and 
8.2  gins  were  used.  The  charge  and  target  were  lowered  to  a  series 
of  depths  and  the  variation  in  amplitude  of  the  2-node  and  3-node 
modes  were  measured  as  the  bubble  period  varied,  due  to  depth,  from 
in  phase  to  out  of  phase  with  the  modes.  The  effect  of  charge 
standoff  was  also  investigated.  For  the  2-node  mode  agreement  was 
again  goou  everywhere  but  for  t.ie  3-noUe  mode  the  experimental 
amplitudes  were  now  greater  than  prodir.  ted.  Prior  to  the  tests 
vibration  measurements  were  made  both  in  air  and  water  tvi  determine 
the  K.odc  sha.es  and  f regencies  and  to  determine  the  values  of  L.. 
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Unlike  the  floating  box  model,  the  non-unifam  muss  distribution  in 

the  cylinder  led  to  minor  differences  in  the  mode  shapes  for  the  two 

cases  so  that,  for  the  cylinder,  condition  (c)  was  violated.  To 

minimise  the  effect  of  vi^ating  condition  (e)  on  these  tests  since 

gravity  migration  could  not  be  balanced  by  surface  repulsion,  the  axis 

of  the  cylinder  was  hung  vertically.  Gravity  migration  was  consequently 

parallel  to  the  axis  and  its  effects  minimised.  In  any  case,  for  such 

tir\y  charges,  particularly  when  fired  at  depth  (between  12  ft  and 

900  ft  in  the  tests),  migration  is  very  small. 

From  the  experiments  it  is  clear  that,  inside  the  restrictions 

of  assumptions  (a)  to  (e),  and  to  some  extent  outside  them,  the  analysis 

is  quite  accurate.  The  good  agreement  also  indie  .tes  that,  at  least 

on  snail  scale,  the  assumptions  (f)  and  (g)  are  justified,  Chertock 

in  fact  give:,  an  argument,  although  not  a  completely  convincing  one, 

why  the  assumption  (f)  should  be  generally  true  even  though  initially, 

when  the  explosion  emits  the  shock  wave,  the  fluid  motion  is  far  from 

incompressible.  The  analysis  however  docs  not  provide  a  fully 

reasonable  method  of  computing  the  whipping  of  ships  generally.  Very 

few  ships  car.  be  regarded  as  proporti  >nal  bodies  and  in  consequence 

the  coupling  coefficients  L.  .  must  be  found  in  some  way.  The 

—  J 

approximate  analysis  only  applies  to  distant  charges  wire  re  the 
distance  fro;  the  charge  to  the  ship  is  large  compared  with  the  ship 
cross  section  and  there  is  no  indication  of  how  rapidly  the  approxima¬ 
tion  breaks  v.own  at  closer  ranges.  For  very  small  charges  of  the 
type  used  in  the  experiment  and  with  the  pressure  still  atmospheric 
or  greater,  explosion  bubbles  tend  to  migrate  very  little.  On  full 


scale  this  ii.  not  true  and  bubbles  migrate  very  vigorously.  Such 
migration  ca:.  lead  to  marked  differcces  in  the  rcsi>onsc  of  ships 
to  charges  vertically  beneath  them  compared  to  charges  at  some  distance 
to  one  side.  In  Chertock's  analysis  conditions  change  only  slowly 
with  stand-off  from  the  ship's  centreline.  Lastly,  being  entirely 
dependent  on  elastic  inodes  of  vibration,  it  is  very  difficult  to 
extend  the  analysis  as  it  stands  to  the  analysis  of  ship  motions 
which  arc  violent  enough  for  plastic  deformation  to  occur. 

These  points  are  all  investigated  in  some  detail  in  the  thesis 
and  a  practical  computation  scheme  is  developed.  Since  methods  of 
calculating  the  normal  mode  frequencies  and  shapes  of  ships  have  been 
well  established  over  the  last  half  century,  their  basis,  strip  theory 
for  the  hydrodynamics  and  an  elastic  free-free  bean  reia-csen tation 
for  the  ship,  is  used  here  to  re-foimulate  the  whipping  equations 
to  anply  to  .ill  ships,  including  those  which  are  not  proportional. 

The  forcing  idion  of  the  bubble  pulsations  on  the  ship  can  then  be 
interpreted  in  terms  of  strip  theory,  which  provides  also  a  veiy  sim.  le 
means  of  including  the  effects  of  bubble  migration.  Use  of  the  basic 
equations  of  an  elastic  beam  also  permits  a  relatively  simple  extension 
to  cover  plastic  hinge  formation  in  the  ship.  The  equations  used 
are  developed  in  the  subsequent  sections. 

To  etuib.c  the  powerful  iiatrix  manipulation  subroutines  readily 
available  in  most  computers  to  be  used  in  computing  the  ship  response 
a  finite  ele. ieiit/l::.pcd  mass  model  has  been  used. 
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C  Th»‘  Elastic  Representation  of  the  Ship 

Since  the  latter  part  of  the  last  century  it  has  been  customary 
to  represent  a  ship  as  an  elastic  beam.  Initially  t'ais  was  largely 
a  necessary  assumption  to  make  strength  calculations  tractable  but 
eventually,  in  1905,  a  series  of  static  bending  moment  experiments .were 
carried  out  on  a  destroyer  !LM5  V/'dl.F  (Cl).  These  confirmed  that  the 
us-.al  thin  c  juations  gave  approximately  the  correct  longitudinal 

stress  distribution  across  each  cross-section,  although  some  deviations 
from  a  linear  stress  distribution  were  apparent  near  the  deck. 

The  results  indicated  the  ’apparent*  modulus  of  elasticity  to  be 
close  to  the  usual'tensile  test'valac  for  low  bending  moments  but 
to  fall  in  value  for  the  higher  moments.  Calculated  deflections, 
when  compared  with  the  measured  deflections  gave  much  lower  apparent 
values  for  ti«  modulus.  The  discrepancy  in  the  deflections  was 
later  shown  to  be  due  to  the  neglect  of  shear  deformations  and  in  1924 
Taylor  presented  an  excellent  paper  (C2)  describing  a  modification  to 
the  straight-forward  beam  theory  which  made  allowance  for  the  effect 
of  shear  lag  on  the  distribution  of  stress  in  beans  of  ship  type 
cross-section.  He  showed  tha'.  under  typical  conditions  shear  could 
make  significant  contributions  to  the  total  shi:  deflection. 

Since  that  tine,  until  very  recently  when  three-dimensional 
finite  clement  stress  analysis  computer  programs  became  available  to 
calculate  stress  distributions  with  some  accuracy,  the  position 
has  not  changed  greatly.  A  number  of  static  bending  experiments 
C3-C8  have  been  carried  out  on  ships  and  all  have  confirm*  tuat  beam 
theory  is  quite  satisfact iry  for  computing  both  stresses  and  dcflecti  ns 
even  to  the  point  of  failure  (C3  and  Cdj  provided  that  long  deck  houses 
are  not  involved.  Such  deckhouses  appear  to  behave  in  a  variety  of 
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ways  depending  on  how  they  are  supported.  Instead  of  the  stresses 
increasing  away  from  the  neutral  axis,  except  in  very  Ions  well 
supported  deckhouses,  the  stresses  decrease.  If  the  deck  on  which 
they  are  supported  is  flexible  enough,  they  nay  even  have  a  curvature 
opposite  to  that  of  the  main  hull.  However,  in  calculating  the 
effect  of  deckhouses  on  longitudinal  deflectiors  it  has  become 
customary  to  treat  the  deckhouse  as  a  normal  art  of  the  beam  (if 
it  is  long  enough)  but  to  assign  to  it  a  reduced  effectiveness  factor 
depending  on  its  size  shape  and  support.  Theories  and  wetnods  for 
calculating  such  effectiveness  factors  arc  given  in  References  C9- 
C',2. 

The  application  of  beai  theory  t<  the  vibration  of  ships  has  a 
history  as  long  as  its  ap;  lication  to  static  calculations  and 
Reference  C ■ 3  gives  an  excellent  review  of  the  work.  So  long  as 
shear  deflections  arc  included  and  an  entrained  mass  of  water  allowed 
for,  excellent  results  can  be  achieved  in  calculating  the  vibration 
frequencies  of  the  first  few  lowest  modes.  However,  shear  deflections 
play  an  increasingly  important  part  for  the  higher  modes  and  the 
distortions  of  plane  sections  (assumed  to  remain  plane  in  normal  beam 
theoiy)  become  important.  Also,  higher  nodes  begin  to  lave  frequencies 
similar  to  some  of  the  transverse  distortion  modes  and  combined  modes 
appear.  i; i:\ally  the  usual  strip  theory  netaod  for  allowing  for  the 
entrained  water  mass  becomes  increasingly  inaccurate  for  the  higher 
modes.  In  cnuscffucncc,  although  high  accuracy  has  been  claimed 
(Cl 4)  for  calculations  of  this  type  up  to  the  10  node  mode,  it  is 
generally  accepted  that  only  the  first  few  modes,  with  up  to  about 
5-mdes,  can  be  calculated  with  confidci.ee.  Fortunately  in  explosion 
induced  whipping  it  is  the  lowest  few  modes  which  turn  out  to  be  the 
most  significant. 
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To  apply  the  lumped  :nass/wcightless  beam  idealisation  to  a 
ship  it  may  be  assumed  divided  into  a  number  of  equal  length  sections. 
If,  as  is  normally  the  case,  the  ship  structure  changes  little  in 
the  length  of  any  section,  it  can  be  assumed  constant  within  each 
section  with  characteristics  corresponding  to  those  at  the  centre 
of  the  section.  The  forces  and  moments  acting  at  the  ends  of  any 
section  are  then  readily  determined  in  terms  of  the  deflections  and 
rotations  there,  all  the  ship  weight,  buoyancy  and  inertial  forces 
oeing  considered  as  lumped  at  the  ends  of  sections.  With  the  sign 
convention  shown  in  Figure  Cl  ,  the  standard  beam  equations  are 


3x 

.2 

M  = 

dx 

If  tlie  angle  of  rotation  at  a  cross-section  is  y»  then  the  slope 

of  the  deflection  curve  is 
dy  S 

<T£  =  Y  T  T7S 

s 

where  A  is  the  area  of  the  cross-section  effective  in  shear,  and 
s 

G  the  shear  modulus.  The  integral  of  these  equations  may  then  be 


expressed  in  the  form 


1  x  m  (£et  ~  rd  [yL 

JL  _x_  { 

El  2EI  Yl 


0  0  1 


0  0  0 


where  the  suffix  L  refers  to  conditions  at  the  left  hand  end  of  the 
beam  section.  If  the  conditions  at  the  right  hand  end  are  substituted 


the  equations  may  be  inverted  to  give  the  forces  and  moments  required 
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at  each  end  cf  the  bean  to  maintain  there  giv~r  displacements  and 
rotations,  Tlu  equations  are 


where 


_  12(1  +v)I  _ 

€  =  ■  ■  ■■  +  1  ■  ■  and  a  = 


•&3(1  +  2e) 


If  there  are  (n  -  1 )  beam  sections  and  the  displacements  (y^*****^) 

and  rotations  (y1»...#yn),  are  imposed  at  the  n  beam  section  ends, 

the  forces  (5  - ...,S  )  (considered  positive  upwards)  and  moments 
l  n 

(M  (considered  positive  anticlockwise)  necessary  to  maintain 

the  imposed  shape  will  be  given  by 


si  fA  3“]  ry 
v  "  /  cj  X 


where  the  elements  of  the  (nxn)  matrices  A,  B  and  C  are  given  by 

2 

a..  ,  s  6a  ;  b..  =  -3».  .1.  .  ;  c.  =  a.  I  (i-e.  .) 

n-1  i-1  ’  n-1  l-l  i-I  n-1  l-l  l-i  x-l 

a..  =  6(a  ,  +  a.)  ;  b..  r  -3a  <  +  3a 

ii  i-1  i  n  i-I  i-I  i  i 

c..  =  a.  ,4?  ,( 2  +  e  )  +  a.4?(2  +  e.) 


'ii  i-1  i-I 


i-l  ii' 


a..  ,  =  -6a.  ;  b..  =  3a.£.;  c..  =  a.£.(i  -  e.) 

ii+l  i  ii+i  i  i  n+1  ii  i 

Here  (e.  a.)  arc  the  values  given  by  (C3)  for  the  beam  section 

and  a  and  c  arc  defined  to  be  zero* 
o  :» 


The  stiffness  matrix  in  equation  (C4-)  may  be  readily  shorn 
to  allow  the  two  degrees  of  freedom  in  the  vertical  plane  (translation 
and  rotation)  expected  of  a  free-free  beam,  without  giving  rise  to 
forces  or  moments.  It  is  therefore  of  rank  (n-2). 

Under  conditions  where  the  applied  moments  are  known  to  be  very 
small  it  is  not  possible  to  specify  the  rotations.  These  will  ther. 
te  related  to  the  displacements  by  the  equation 

Y  =“C",BTy  (C5) 

and  the  forces  necessary  to  maintain  the  displacements  y  are  then 
£  =  [A  -  wfVly  (C6) 

This  form  can  be  useful  since  for  ship  vibrations  the  rotary  inertias 
are  generally  small  enough  to  be  neglected. 
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D  The  'lydrudvn.  mic  Farces  on  the  Ship 

Since  'strip'  theory  lias  over  the  l  ist  x  years  been  the  principal, 
in  fact  almost  sole  method  of  allowing  for  the  entrained  mass  of  water 
for  vibration  problems  in  still  water,  it  has  been  used  here  in  the 
initial  formal. ition  of  the  equations.  An  alternative,  more  accurate 
method  is  presented  later.  Strip  theory  was  first  introduced  in 
1929  by  Lewis  (D » ;  and  has  changed  very  little  since  that  time. 

Lewis  exploited  the  long  slender  nature  of  ships  by  assuming  that 
close  to  the  ship  tire  liydrcdynamic  flow  would  be  largely  confined  t<» 
planes  perpendicular  to  the  ship  axis  and  so  at  each  cross-section 
would  approximate  closely  to  the  two-dimensional  flow  around  a  cylinder 
having  the  same  cross-section  as  the  ship  hull  at  the  point  of 
interest.  This  should  clearly  give  a  good  approximation  to  the 
local  flow  except  near  nodes,  where  the  ship  motion  is  largely  a 
rotation  about  the  node,  or  where  the  ship  cross-section  is  changing 
fairly  rapidly,  for  example  near  the  bow  and  stern.  Fear  these 
exceptional  points  there  will  be  a  significant  flow  along  the  ship 
axis.  Ignoring  these  points,  with  their  associated  three-dimensional 
flow  fields,  the  force  per  unit  length  on  any  cross-section  of  the 
hull  will  be  identical  to  that  of  the  corresponding  two  dimensional 
flow  on  the  cylinder.  Conformal  transformations  were  uspd  to 
extend  the  well-known  flow  around  a  two-dimensional  circular  cylinder 
into  flows  around  appropriately  shaped  cylinders.  Lewis  use’  the 

Q  k 

transformation  Z  =  a  +  -  + 

8 

to  give  the  flow  around  the  cross-sections  shown  in  Figure  Dl , 

The  added  masses  per  unit  length  (p)  were  related  to  the  added  mass  of 
a  circular  cylinder  of  the  same'beam'as  the  ship  section  by  the  formula 

p  =  C 
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The  coefficients  C  are  given  for  the  cross-sections  in  Figure  1. 

In  order  to  allow,  at  least  partially,  for  the  three-dimensional 
flow  effects,  Lewis  compared  the  strip  theory  flow  around  a  prolate 
spheroid,  vibrating  in  a  certain  'mode'  shape  similar  to  a  ship 
vibration  Mode,  to  the  exact  flow  deduced  by  an  expansion  of  the 
flow  m  terms  of  spheroidal  harmonic  functions.  The  kinetic 
energy  of  the  approximate  flow  is  naturally  greater  than  that 

of  the  exact  flow  (Kelvin's  minimum  energy  theorem) 

so  tliat  the  deduced  added  mass  will  also  be  greater.  Lewis 
suggested  that  the  calculated  added  mass  distribution  for  the  ship 
shape  should  be  reduced  by  a  factor  J,  the  ratio  of  the  exact  to 
the  a  iprcximate  kinetic  energies  for  a  spheroid  with  the  sane 
beam  to  length  ratio  as  the  ship. 

The  whole  procedure  g.vcs  very  good  results  for  the  lower 
frequency  vioration  nodes  of  ships  and  forms  the  basis  of  almost 
all  current  ..etheds  for  calculating  the  effect  of  the  water. 

It  has  two  rather  basic  disadvantages  however.  Firstly,  the 
reduction  factor  J  for  the  3-D  flow  effect  depends  on  the  shape 
of  the  node  :c  be  calculated,  which  ra  :cs  it  difficult  to  formulate 
a  single  o«,u  ition  whose  eigenvalues  art  the  ship  frequencies; 
u  different  equation  has  to  be  solved  for  each  mode.  Secondly 
there  is  no  guarantee  the-  .n  exact  reduction  factor  would  not 
vary  along  the  length  of  the  ship,  tr!  in  general  this  is  the  case. 

Although  the  method  is  still  the  basis  for  entrained  mass 
calculations,  numerous  minor  improvements  have  been  :uuc  since  Lewi:;' 
paper.  'fay  .or  (D2),  arrived  at  the  sane  a  proximate  strip  theory 
independently  of  Lewis  and  again  calculated  i  red -c'. inn  factor  by 
comparison  w„th  the  solution  for  a  spheroid.  now  ever,  fv/lor's 
solution  for  the  spherical,  vibration  allowed  a  bending  type  of 


motion,  as  well  as  the  shearing  motions,  on  which  Lewie  calculations 
were  based,  in  the  boundary  condition  on  the  spheroid.  In  consequence 
his  reduction  factors  give  smaller  entrained  nesses  than  those  of 
Lewis.  Curves  for  both  types  of  reduction  coefficient  are  shown, 
for  the  2-node  vibration  mode,  in  Figure  D2.  Taylor  also  gave  a 
second  so:  of  reduction  coefficients  based  on  an  exact  solution 
for  the  three  dimensional  motion  of  an  infinite  cylinder  vibrating 
transversely,  its  velocity  distribution  varying  sinusoidally 
along  its  length.  These  reduction  factors  arc  also  included  in 
Figure  D2.  To  calculate  them  the  Tnylor  Form  of  ellipsoid  motion 

i. 

was  assumed.  This  lias  nodes  at  hk  L(1-23Vl  )  /2/5»  measured  from 

its  centre.  The  distance  between  these  nodes  was  assumed  to  be 

the  half  wavelength  for  the  cylinder  motion  whence  the  appropriate 
reduction  factor  can  be  readily  found.  That  they  agree  slightly 
better  with  Lewis*  original  results  is  not  surprising  since  the 
boundary  conditions  used  by  the  two  methods  coincide  when  the 
vibrating  body  has  a  constant  cross-section.  Only  when  appreciable 
changes  in  section  occur,  as  close  to  ends  of  the  spheroid-  or,  to 
a  lesser  extent,  near  the  bow  and  stern  of  ships,  do  the  boundary 
conditions  and  hence  the  calculated  reduction  factors  differ. 

Froimska (D3), considerably  generalised  Lewis'  two-dimensional 
cross-sections  by  considering  the  truns formation 

Z  =  »  +  a/ftm  +  b/^n. 

for  the  enscs  (n,n)  =  (1,5),  (1,7)  and  (3,7).  Lewis*  sections  of 
course  are  for  (m,  n)  =  (l,  3).  From  these  results  Prohaska  found 
that  for  sections  with  some  parts  concave  outwards,  the  added  mass 
was  very  close  to  tha'  for  the  same  section  with  the  concavity 
replaced  by  a  tangent  line.  taking  use  of  this  device  he  then 


found  the  coefficient  C  could  be  determined  with  good  accuracy  at  a 
function  simply  of  t lie  section  area  coefficient^  (the  underwater  area 
divided  by  b;an  x  draft,  BD)  and  the  ratio  beany^draft.  He  gave 
the  functional  form  as  the  graph  reproduced  in  Figure  D3.  This 
graph  enables  the  added  mass  to  be  found  to  reasonable  accuracy 
without  the  tedious  comparison  of  the  actual  cross-section  with 
those  of  Lewis  to  find  the  best  fit. 

The  set  of  cross-sectional  shapes  fo.r  which  the  added  mass  is 
known  has  been  fur  the;  extended  in  recent  years  by  Landweber  and 
Macagno  and  co-workers  1)4,  D5.  These  authors  used  a  three  parameter 
fit  for  the  sections  rather  than  the  two  para»,ietcr  fits  used  by 
Lewis  and  Prohaskn,  Lacagno  also  made  use  of  ’closc-fit*  conforn*l 
mappings  of  sections,  and  other  authors,  e.g.  Hof f man (DC) have 
demonstrated  similar  methods.  These  use  an  arbitrary  number  of 

coefficients  in  the  series 

_  -\  ®3  a 

Z  =  z+-i+-*  +  -5.  +  ••• 

*  z3  »J 

to  obtain  ns  close  a  fit  to  the  actual  section  as  desired.  Quite 
complicated  computer  programs  have  to  be  used  however  to  obtain  the 
coefficients  for  the  best  fit  transformation. 

.Vendel  (r?)  has  given  the  added  mass  for  a  squar-'  cross-section 

f) 

with  •bilge-icecls’ ;  plates  projecting  at  45  from  each  comer. 

These  results  are  particularly  interesting  as  they  appear  to  be  the 
oniy  ones  available  on  the  effects  of  such  keels.  The  method  used 
to  obtain  thorn  was  on  application  of  the  Schwarz-Christof fel  transforna* 
tion.  Unfortunately,  the  method  is  very  tedious  to  apnly  nnd  the 
results  were  not  extended  to  general  rectangles  with  bilge  keels. 

The  results  for  such  shapes  can  however  probably  be  estimated 
fairly  well  by  Prohaska’s  method;  simply  using  the  added  mass  for 
the  recta agio  which  circumscribes  the  ends  of  the  bilge  keels 
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(the  dotted  rectangle  in  Figure  D4).  For  the  square  section  fc'endel 
quoted  the  results  in  the  table  for  the  increase  in  the  coefficient 
C  for  various  values  of  d/b  (the  bilge  keel  width  to  half  bean  ratio). 
If  the  added  nass  of  the  circumscribed  section  is  taken,  the  factor 
by  which  C  is  increased  (relative  to  the  original  rectangle)  is 

*  2 

(i  +  *  The  ■'  increases  due  tc  this  factor  have  been  included 

in  the  table. 


Wendel's  Bilge  Keel  Effect 


The  approximate  result  should  be  close  enough  for  practical  purposes. 
V/endel  also  commented  on  the  constancy  of  the  3-D  flow  correction 
factor  along  the  length  of  the  ship  and  showed  that  in  fact,  even 
for  such  an  unlikely  case  as  a  rigid  sphere,  the  factor  can  be 
exactly  constant,  at  least  for  rigid  body  motions. 

This  point  was  also  considered  by  Kaplan(DS) who  shoved  that  for 
arbitrary  'slender'  bodies  the  reduction  factor  is  constant  at 
least  for  rigid  body  motions.  He  also  showed  that  the  reduction 
factor  is  constant  for  both  Ltuis'  and  Taylor’s  assumed  2-node 
vibration  modes  of  a  spheroid.  The  proof  for  the  Taylor  mode  shape 
however  Is  not  quite  correct  and  in  fact  the  local  reduction  factor 
does  vary  for  this  case.  Kaplan’s  paper  gives  an  excellent  survey 
of  tie  virtual  mass  problem  cs  it  stood  In  1959*  In  it  he  discusses 


t  ic  effects  of  viscosity,  compressibility  and  surface  wes  and 
concludes  th..t  all  sliould  be  very  scrJl  (except  possibly  viscosity) 
for  normal  vibration  frequencies. 

The  surface  wave  problem  has  been  considered  by  a  number  of 
authors,  e.g.  Ursell(DS),  Porter(D. 0)nnd  Choung  ,*!oo<:  Lee, (Oil). 

The  results  all  show  that  for  the  pitching  and  heaving  ship  nodes 
the  effect  of  the  surface  waves  changes  the  added  mass  appreciably 
and  provides  significant  damping.  However,  for  the  vibration  Bodes, 
Which  are  of  rather  higher  frequency,  the  effect  of  surface  waves 
becomes  negligible. 

Hydrodynamic  Forces  in  Moving  water. 

According  to  strip  theory  the  hydrodynamic  force  on  a  ship 
cross-scctin'i  which  has  an  acceleration  v  is 


-fiv  per  unit  length  where  p  =  CJ 


(01) 


B  is  the  ship  beam,  p  the  density  of  water,  C  a  factor  derived  from 
a  2-diiaensional  flow  problem  and  J  is  a  correction  factor  for  3-D 
flow  effects.  To  help  detemnine  the  force  when  the  surrounding 
water  is  also  in  motioi.  it  is  helpful  to  shor  the  dependence  of  C 
on  the  potential  of  the  2-D  flow  problem.  With  the  geometry  of 
Figurep5  the  potential  <f>  satisfies  the  following  equations:- 

-  o  (incompressible  flow) 

<t>  =  0(b),  m«o  (only  a  dipole  field  at  large  distances) 

££  =  -  v  cos 0"*  on  T  (the  shio-shnned  cylinder  cross-section) 
on 

If  the  potential  for  Hr*  case  v  =  j  is  denoted  by  <p,  then  #  =  v  # 
and  <f>  satisfies  y2^  _  # 

♦  =  #(  -) »  r  ■»» 


2/, 


-  cjs6  on  1’ 


flic  usual  fora  of  Dcrnnuilli's  c -uatiul.  gives  the  •  rrssurc  in  the 


fluid 


l3  =  I30  +?lt(v^)  "  ;^U 


;,1  cr .  :  is  t’ac  IrVrostatic  pressure  ;jiu  u  LI :c  lo':  J  velocity# 
o 

h'cglcctinj  pu“  compared  to  p  ^  jives  ;>  =  7,  +  pv^i  .-.mi  the  force 

actinj  on  the  cylinder  is  :>i..v>ly 


-jr  d»  -  -pvj  <f>dn 


so  that,  co.  ..>a ring  v.ith  (Dt), 


Jr 


The  factor  _  allows  for  the  u  per  half  of  the  fluid  not  being 
present  since  y  =  0  is  the  free  £>•  face. 

..lien  the  surrounding  fluid  also  has  1  velocity,  U ,  ut  infinity, 
the  potential  f>'  must  satisfy  the  equations 


V“^  =  o 

if>'  =  -Uy  +  0(p),  >*  -* 


N  I  *  * 

-r*  -  -v  cos6  ?r\  I 

on 

lleariy  this  potci.ti.il  is  given  by  <f>'  -  -Uy  -  (U-v)<£  for  this  s-  tisfici 

the  first  tv.o  conditions  and  on  T 

^  =  -  Ur^“  -  (  f-v)^  =  -  V  cosO  *  -  (L-v)  ( ~  cosO*)=  -v  cos  0 
dn  On  dr  ' 

as  repaired. 

The  pressure  it.  t.<c  •.  at.r  is  no.,  (again  neglecting  pu“  L.-r.ss) 


13  a  :  o  *  pit c"  !jy  "  (u  "  v)^3 


and  the  fore  •  on  tlx-  cylinder  is 


-4  pd»  =  pU^yu#  ♦  p(U  -  v)4y  tli 
j  v  J  j 
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Nov  £  yd®  is  the  ,r.a  within  the  contour  T,  i.e.  2A  -..here  \  is  the 
underwater  area  of  the  cross-section,  so  that  the  force  on  the 
underwater  section  is 

f  =  pAU  ♦  |j.(  U  —  v)  (D6) 

which  reduces  to  the  previous  case  when  U  =  o.  This  equation  is 
in  fact  rather  obvious.  The  second  term  is  the  inertial  force  due 
to  the  relative  acceleration  of  the  cylinder  and  the  fluid  and  the 
first  is  the  'buoyancy'  force  due  to  the  pressure  gradient  required 
ir.  the  water  to  give  it  the  acceleration  U.  Under  gravity  the 
fluid  effectively  has  an  acceleration  l'  =  g  and  the  first  tern 
becomes  the  weight  of  the  ■  ispu.ccd  water,  as  is  assumed  in 
hydrostatic  calculations. 

This  result  is  applied  in  the  sane  way  as  normal  strip  theory. 

The  vertical  velocity  of  the  water  due  to  the  explosion  bubble  notion 
can  be  calculated  at  each  point  of  the  ship  axis,  (the  intersection 
of  its  centreplane  with  the  waterline  plane),  as  though  the  ship  were 
not  present.  If  the  charge  is  far  enough  from  the  ship  for  this 
velocity  distribution  to  vary  only  slowly  along  the  length  then 
at  cacli  cross-section  the  disturbance  potential  necessary  to  correct 
the  local  flow  for  the  presence  of  the  ship  will,  as  above,  be 
approximately  -(U-v)<£  being  appropriate  to  the  section  sha-e) 
and  the  force  or.  the  ship  will  be  as  given  by  (06). 

Since  pi  -os  been  expressed  in  the  fora 


it  is  coiiven.cnt  to  express  Ap  si.ii1.  rly  so  that 
Ap  = 


Id  this  case  however*  since  the  ’buoyancy'  force  depends  on  a 
pressure  gradient  independent  of  the  shape  of  the  ship,  no  three- 
dimensional  correction  factor  is  retired.  The  values  of  the 
coefficient  C  for  Lewis'  sections  have  >>cen  added  to  Figure  D1 . 

C  is  related  to  the  usual  section  area  coefficient,  p,  by 
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VALUE  OF  INERTIA  COEFFICIENTS  C  AND  C 
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THREE  DIMENSIONAL  REDUCTION  FACTOR 


FIGURE  D2 
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WENOEL'S  RECTANGLE  WITH  BILGE  KEELS 


FIGURE  D4 


GEOMETRY  FOR  ACCELERATING  FLOW 
PAST  A  2  -  D  SHIP  SECTION 


FIGURE  D5 


33 


E  Equations  of  Motion 

The  only  external  forces  acting  on  the  ship  are  the  hydrodynamic  ones 
discussed  in  the  previous  section.  These  form  a  continuous  distribution 
of  force  along  the  ship.  In  addition  to  these,  the  inertial  forces  due  to 
the  ship's  own  mass  also  act  and  again  form  a  distribution  along  the  entire 
length.  The  fonr—lation  of  the  elastic  nature  of  the  ship  assumed  only 
forces  acting  at  tr.e  discrete  'nodes'  between  the  beam  sections.  Current 
methods  of  finite  element  analysis  generally  use  consistent -mass  matrices 
derived  from  energy  considerations.  These  divide  the  mass  between  the 
element  nodes  in  t;<e  most  effective  manner.  However,  for  the  beam 
displacement  shapes  inherent  in  the  elastic  formulation  the  split  is  rather 
complicated  and  al^o  much  of  the  ship  mass  comes  from  material  such  as 
machinery,  transverse  bulkheads,  etc.,  which  is  not  continuously  distributed. 

It  is  consequently  convenient  and  reasoiable  to  simply  split  the  ship  into 
n  equal  sections,  lump  the  total  mass  and  inertia  of  each  section  into  a  single 
mass  and  inertia  at  the  centre  of  the  section,  and  take  the  elastic  beams 
to  be  the  (n-l)  equal  lengths  connecting  these  masses.  The  length  t 
of  each  elemental  oeam  is  L/n.  If  the  variation  in  cross-sectional  shape 
is  fairly  slow,  the  lumped  hydrodynamic  force  and  moment  at  the  ith  mass 
will  be  from  (D6)  (pi  again  being  the  added  mass  per  unit  length). 


shdi  *  ['*i(Vyi)  *  Vi1* 

•  .*•  n  * 

u  r^i,:  ",  ,d“i  dt*i  •  -  dUil  <3 

Mhdi  dx  ui  yi  +  dx  ~  dx  +  dx  ui  +  dx  12 

2 

Ci  d 

where  all  terras  of  order  — 2  and  ( j-)  have  been  neglected. 

dx  X 

The  corresponding  hydrostatic  terms  (loss  of  buoyancy  due  to  upward 
displacement  will  be 
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The  normal  ii.rtial  force  and  moment  .ire 


S_  .  -  -m.y . 
Ii 


MIi  =  “ri  dx 


Combination  of  ions  (J-l  !  to  (C3)  gives 


S.  =  -(a.+ra  }y.  -  k.y.  +  (m  .+5  )u 
i  l  wi  i  ii  wi  wi  i 


M .  =  -r  y.  -  k' .  y  . 
l  wi  l  ri  l 


(r'.+r'.)u.  -  ( r .  +  )y . -k  Y .  + 
tvi  wi  i  i  wi' 1 1  r1  i 


( r  .  +  r  . ) — t— 
wi  wi  dx 


•here  m,.  =  ».l,  ^  =  p.C.  ^  =■  ft  "„/•  ?„i  =  k  "wi*2 


ki  =  V<  kri  * 


r*  n  -L  e3  r  -  L  e: 

rwi  1 2  dx  ^  *  rwi  "  ) 2  dx' . 


i  dL-  ^ 

k'  _  J _ i  P 

ri  "  12  dx  ^ 


In  matrix  form  th<  -e  equations  (E4)  become 


M  +  M 


°1  ry  k  o  ~|  r«w+i  0  u"j 

—  ^  W  W  ^  , 

••  *  ’ 

R+RW  r  "  k'  K  y  +  r'+r'  R  +r  u'  i 

1 _ L  I  [_  j  j_  ]_'j 


;.'i  cm) 


where 


U  =  ox  and  tl,e  s,lbi,iatrices  arc  all  diagonal  with  the  obvious 


elements.  This  equation  is  slightly  inconvenient  as  it  involves  the 
off-diagonal  submatrices  R',  IT  and  k'.  Over  most  of  the  length  of  a 
ship  the  cross-section  is  almost  constant  so  that  these  terms  will  be 
very  small  except  possibly  near  the  ends.  Strip  theory  already  ignores 
longitudinal  flow  due  to  changing  cross-section  in  these  regions.  It  is 
therefore  consistent  with  strip  theory  to  neglect  the  terms  and  use  the 


equation 


On  applying  O' Alembert  s  Principle,  these  arc  the  forces  and  moments  neccs^-iry 
to  maintain  the  displacement  shape  (y,  y)  in  the  elastic  equations  (C4)» 

•v  «v 

The  equations  of  forced  motion  of  the  ship  are  therefore: - 


and  for  the  still-water  vibration  the  right  hand  side  may  be  dropped. 

In  general,  for  ship  vibration,  it  has  been  found  that  rotary  inertia 
effects  are  usually  fairly  small.  On  neglecting  the  appropriate  terms  and  putting 


and 


(M  ♦  M  )y  +  (A+K  -  BC~lBT)y  =  (H  +M  ,)u  (E8) 

W  -V  -X-  W  W  ~ 

Neglect  of  rotary  inertia  in  these  lumped  mass  equations  of  motion  differs 
from  its  neglect  from  the  normal  Timoshenko  beam  equations  since  here  it  has 
a  slightly  different  meaning.  In  the  Timoshenko  beam  equation  rotary 
inertia  refers  to  the  inertia  of  an  infinitesimally  thin  cross-sectional 
slice  of  the  beam,  and  so  is  related  purely  to  the  depth  of  the  beam  and  the 
amount  of  material  away  from  the  centre  of  rotation.  Here  the  rotary  inertia 
includes  this  term  but  also  includes  a  term  due  to  part  of  each  lumped  mass 
being  a  finite  distance  from  the  mass  in  the  longitudinal  direction.  It 
therefore  depends  on  the  distance  between  the  lumped  masses  and  so  is  a 


function  of  the  beaa  representation  rather  than  an  intrinsic  property 
of  cross-sections  of  the  beaa.  Only  if  the  nuaber  of  luaped  masses 
representing  the  beaa  were  allowed  to  tend  to  infinity  would  the 
two  foras  of  rotary  inertia  coincide.  In  the  equations  used  here, 
each  fluid  rotary  inertia  is  entirely  due  to  this  limping  of 
longitudinal  forces  into  point  load  and  so  vanishes  entirely  in 
the  liait  of  an  infinite  nuaber  of  lumped  aasses. 

It  aay  be  noticed  that  if  the  ship  is  in  still  water  and  is  not 
vibrating  then  the  equations  (E7)  or  (£8)  are  satisfied  by  zero 
deflections  everywhere.  They  therefore  do  not  include  any  static 
bending  aoaents  and  shear  forces  due  to  mismatches  between  the 
hydrostatic  buoyancy  and  weight  distributions.  These  static  terms 
were  oaitted  in  equation  (E4)  for  siaplicity  in  the  dynamic  equations. 
Ai\y  stresses  deduced  froa  equations  (E7)  and  (E8)  during  aotion  are 
to  be  superimposed  on  the  existing  static  stress  distribution  in  the 
ship. 

The  equations  include  the  buoyancy  terms  K  and  and  so  contain 
the  heaving  and  pitching  ship  frequencies,  nt  least  for  small 
displacements  where  the  additional  hydrostatic  buoyancy  forces  are 
linear  functions  of  the  displacement* 
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The  Normal  Mode  Equations 

The  equations  (E7)  and  (E8)  are  both  linear  and  of  the  general 

fora 

HJ  t  Sy  =  M  (Fl) 

1  -v 


where  and  are  diagonal  matrices  and  S  is  a  symmetric  positive 
definite,  stiffness  matrix.  For  these  equations  to  predict  the 
transverse  snip  motions  with  any  accuracy,  the  ship  has  to  be  split 
into  at  least  20  lumped  masses  and  the  direct  integration  of  the 
equations  becomes  rather  tedious.  This  is  especially  so  since  the 
equations  will  then  contain  at  least  20  normal  mode  frequencies  the 
highest  of  which  will  be  little  less  than  that  of  the  vibration  of 
a  single  lumped  mass  on  its  two  weightless  beams,  all  the  other 
masses  being  considered  stationary.  Clearly  such  a  mode  will  have 
no  physical  meaning,  but  its  frequency  will  be  very  high,  and  will 
become  higher  as  the  number  of  masses  is  increased.  Direct 
numerical  integration  of  equations  (Fl)  requires a  time  increment 
no  greater  than  about  l/6th  of  the  period  of  this  highest  mode 
if  numerical  instabilities  are  to  be  avoided.  In  consequence  the 
direct  integration  involves  the  integration  of  a  large  number  of 
equations  using  a  very  small  time  increment. 

A  simple  alternative  is  to  split  the  equation  into  its  normal 
modes  and  to  reject  at  the  outset  the  higher  modes  known  to  have  no 
physical  basis.  Since  is  diagonal  the  simple  transformation 


1/  j/ 

(where  M  ^  is  diagonal  with  elements m/2)  reduces  (Fl)  to 


1,  I- 

eu"*2 


where  E  =  SI!, 


Since  S  is  symaetric  and  positive  definite,  E  will  have 

similar  properties  and  its  eigenvalues  will  be  real  and  positive. 

The  eigenvectors  and  eigenvalues  of  E  can  therefore  be  written 
2 

in  the  fora  it.,  u.  and  the  vectors  will  fora  an  orthogonal  set, 

i.e.  »T  a.  =  M.6 .  .  (F3) 

^  i  »J 

where  li.  is  a  positive  non-zero  scalar.  The  vectors  a.  Fora  a 

V  • 

complete  set  so  that  the  vectors  •  and  Kp  ll2u  ■*>  expressed  as 


jN  /  \  N 

»  =  ,  ,Vl,,i  ;  U  li -U  szp.(t)i., 

~  i=1  ~  1  2~  .  ,  x 


when  E  is  of  order  (N  x  N) * 

Substitution  of  these  expressions  into  equation  (F2)  and  use 
of  the  orthogonality  of  the  a.  gives 

m  2 

ai  *  =  =  •••» 

These  equations  are  independent  and  only  the  few  low  frequency 
modes  of  interest  need  be  solved. 

The  ship  deflection  shapes  corresponding  to  the  normal  modes 
are  given  by 

y.  -  C*  .. 

Li  i 

In  terms  of  these  displacement  shapes  the  orthogonality  relations 
(F2)  become 


li  *iij  =  “isij 
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and  a-U)  +  ^a^t)  =  ~  2  ^P.tfki^k  » =  X*  (f6> 

which  gives  the  equations  for  the  ith  node  coefficients  a^(t) 

in  terns  of  the  hydrodynamic  mass  distributions,  the  displacement 
shape  of  the  ith  node  and  the  fluid  accelerations  at  all  the 
lumped  mass  positions* 

If  the  coefficients  a^{t)  are  found  by  integration  of 

equations  (Fo)  the  actual  ship  displacements,  stresses  etc.  are 
easily  determined  from  the  equations: 

N 

displacement  y(x,t)  =  |_j<i  (Oy^x) 

N 

velocity  £(x,t)  =  2  &.(t)y.(x) 

i=l  1  1 


N 

bending  moment  M(x,t)  =  2  a.(t)M.(x)  (F7) 

i=l  x 


deck  stress  a(x,t)  =  «(x,t)  /  y  =  distance  of  the 

deck  above  the 
neutral  axis  ) 


N 

shear  force  S(x,t)  =  2  a.(t)S.(x) 

i=1  1  1 
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where  y^(x),  M^(x)  and  S^(x)  are  the  displacement,  moment  and  shear 

force  at  the  position  x  when  the  ship  is  in  the  ith  mode  (unit 

amplitude).  With  the  normal  mode  shapes  known  (y..,  y  .  known  at 

ki 

the  kth  mass)  the  displacement,  moment  and  shear  force  for  a  point 
at  a  distance  x  =  x-x^  from  the  left  hand  end  are  given  by 


y.U) 

1  X 

-2 

X 

1 

X 

K 

* — ' 

_ 1 

2EI. 

k 

6EIk 

yki 

«i(x) 

= 

0  0 

1 

-X 

“l 

X 

•H 

C/5 

0  0 

0 

1 

i 

j 

tn 

_ 1 

and  S 

L 

are  given  by 

equation  (C2). 

(F8) 
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Use  ot*  the  Modal  Equations  for  Free  Vibration 


To  illustrate  the  use  of  the  modal  equations  for  free  vibrations 
and  to  estimate  their  accuracy  when  used  in  ship  calculations,  three 
examples  are  considered: 

a.  A  uniform  solid  beam  in  air  and  water 

b.  A  hollow  rectangular  box-beam  in  air  and  water 

c.  A  full  scale  ship 

In  some  of  the  previous  vibration  calculations  in  the  literati  re  ships  havebem 
divided  into  up  to  40  sections  but  20  is  the  more  usual  number  and 
has  usually  been  found  acceptable  for  the  low  modes.  The  examples 
are  therefore  based  on  20  lumped  masses.  The  calculations  were 
carried  out  using  the  computer  program  I  described  later  in  section 
I. 


a.  Uniform  Solid  Beam 

In  1928  in  an  early  series  of  experiments  on  the  effect  of 

entrained  water  on  the  vibration  of  ships,  Moullin  and  Browne 

(G1 )  used  a  series  of  uniform  rectangular  mild  steel  bars.  The 

bars  were  78  inches  long,  2  inches  deep  and  of  varying  thickness. 

The  principal  thickness  was  1  inch.  The  cross-sectional  moment 

of  inertia  for  such  a  beam  is  -BD  (see  Figure  G1 )  i.e.  l/6  in\ 

The  shear  area  for  a  rectangular  section  is  2/3rds  of  the  cross- 

4.  2 

sectional  area,  i.e.  .  Then,  for  representing  the  beam  by 
20  equal  lengths  the  equations  of  section  C  give:- 


l  =  3.9";  I  =  .1667  in4;  Ag  =  1.333  in2 
E  a  30.10**  psi;  v  =  0.3 


€ 


12(Uv)I 


s 


2EI 


r (i  +  2c) 


=  0.1342.1  Oblt/in. 


=  0.1282 


a  = 


The  matrices  A,  B  and  C  are  then  the  (20  x  20)  matrices 


A  =  10 


.805  -.805  0  0  ...  0 

-.805  1.610  -.805  0  ...  o 

0  -.805  1.610  -.805  0  ...  0 


lb/in 


-.805  1.610  -.805 


C  =  \<f 


0 

0  -.805 

.80 

.570 

1.570 

0 

0 

0 

■l  .570 

0 

1.570 

0 

0 

0 

-1 .570 

0 

1 .570  .. . 

0 

• 

• 

• 

• 

c 

• 

-1 .570  0 

1.: 

0 

0  -.1570 

-1 .! 

4.35 

1.781 

0 

0  ... 

0 

1.781 

8.70 

1.781 

0 

0 

0 

1.781 

8.70 

1 .781 ... 

0 

• 

• 

• 

• 

• 

• 

1 .781  8.70 

1.71 

0 

0  1 .781 

4.3! 

theso 

define  the 

shear 

forces  S.(lb)  and  benoin 

lb-in 


M^(lb.in)  for  a  given  displacement  shape  (in).  y.  (radians) 
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according  to 


Since  the  density  of  steel  is  .284  lb/in^,  the  mass  of  each 
'lumped  mass*  is  2.2l6  lb.  For  the  rectangular  cross  section  one 
of  the  curves  in  Lewis'  paper  (G2)  gives  C  =  1 .36! .  The  uncorrected 
added  mass  for  each  section  is  therefore 

b2 

*P(|  )ec  =  o.6oi  ib 

(t.he  usual  factor  ^  has  been  dropped  since  for  a  fully  submerged 
bar  no  surface  correction  factor  is  necessary).  For  a  beam  with 
L/B  =  39*  the  3-D  reduction  factor  is  0,983  so  that  3-D  flow  affects 
the  bar  very  little.  The  buoyancy  term  is  not  actually  needed  for 
the  free-vibratj.on  calculation  but  for  this  example  it  is 
pA4  =  2oBDC  =  .2815  ib.  The  immersion  force  ^  is  of  course  zero  for  a 

submerged  bar.  The  rotary  inertia  of  each  section  about  its  centre 

-g-  BD4(£2+4D2)  =  2.99  lb. in2.  The  rotary  inertia  of  the 
2 

surrounding  water  is  l/l2m and  of  the  buoyancy  it  is 
l/t2mw*  so  that 

Bi  =  2.216  lb;  mwi  =  .591  lb;  m^  ^  .282  lb 

r.  =  2,99  lb. in.2;  r  .  =  .749  lb. in.2;  r  .  =  .358  Ib.in.2 

X  WI  wx 

Using  these  values,  the  normal  mode  analysis  described  in  section 
F  is  readily  carried  out  using  standard  computer  matrix  manipulation 
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subroutines 


For  the  beam  vibrating  in  air  alone  (a  .  -  r  .  =  0), 

wi  wi 

and  with  shear  deflections  suppressed,  (As  made  very  large)  the 
calculated  frequencies  are  compared  with  the  exact  frequencies 
for  a  free-free  uniform  beam  in  the  following  table.  The  exact 
frequencies  are  given  by 

where  x.  has  the  values 
1 

4.730,  7.353,  10.995,  14.137,  17.278,  20.420,  23.561,  26.703, 
29.845,  32.986,  ...  (the  roots  of  cos  x  cosh  x  =  l) 
and  these  frequencies  do  not  include  the  effect  of  shear  deflections 

COMPARISON  OF  EXACT  AND  LUMPED  MASS  FREQUENCIES 
FOR  A  UNIFORM  BEAM  VIBRATING  IN  AIR 


Number  of  nodes 

2 

3 

4 

5  6 

Exact  frequency 

34.13 

94.1 

184.4 

304.9  455.5 

Calculated  frequency 

34.08 

93.6 

182.4 

299.0  442.2 

%  error 

0.15 

0.53 

1.08 

1 .94  2.92 

Number  of  nodes 

7 

8 

9 

10 

Exact  frequency 

636.2 

846  o9 

1088 

1359 

Calculated  frequency 

610.1 

801 .3 

1013 

1243 

%  error 

4.10 

5.39 

6.90 

8.54 

The  agreement  to  within  0.1 5?S  for  the  basic  two  node  frequency 
is  quite  reasonable  and  confirms  the  general  accuracy  of  the  lumped- 
mass  formulation  and  the  computer  routines  used.  For  the  higher 


modes  the  accuracy  naturally  falls  off  but  that  it  should  be  only 
8 in  error  for  the  ten  node  frequency  is  surprising  considering 
that  only  a  20  lumped  mass  system  was  used  so  that  there  are  only 
two  lumped  masses  per  half  wavelength. 

To  show  the  effects  of  rotary  inertia,  shear  deflections  and 
immersion  in  water  on  this  beam  a  series  of  calculations  were  run 
to  produce  the  results  in  the  table.  The  'exact*  results  again 
contain  no  correction  for  shear  effects. 

FREQUENCIES  FOR  A  UNIFORM  BEAM  IN  AIR  AND  WATER 


Number  of 
^''v^nodes 

Case 

2 

3 

4 

5 

6 

1A 

34-1 

94.1 

184 

305 

456 

2A 

34-3 

94.7 

186 

308 

462 

3A 

34.1 

93.6 

182 

299 

442 

4A 

34.3 

94.6 

186 

307 

458 

5A 

34.1 

93.5 

182 

298 

439 

6a 

33.5 

- 

- 

- 

- 

1W 

30.3 

83.6 

164 

271 

405 

2W 

30.4 

84.1 

165 

274 

410 

3W 

30.3 

83.2 

162 

266 

393 

4* 

30.4 

84.1 

165 

273 

407 

5W 

30.3 

83.1 

162 

265 

391 

6W 

30.2 

*82 

- 

- 

- 

Cases:  l  ;  'Exact'  uniform  beam  frequencies  (without  shear) 


2  ;  Calculated  excluding  shear  and  rotary  inertia 


3  ;  Calculated  excluding  shear  but  including  rotary  inertia 

4  ;  Calculated  including  shear  but  excluding  rotary  inertia 

5  ;  Calculated  including  both  shear  and  rotary  inertia 

6  ;  Measured  value 
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Clearly  neither  rotary  inertia  nor  shear  affects  the  frequencies 
greatly  for  this  beam.  Shear  in  fact  has  almost  no  effect  but 
inclusion  of  rotary  inertias  drops  the  6-node  frequency  by  about 
4£o.  The  agreement  with  the  measured  frequencies  is  within  the 
experimental  error. 

The  main  conclusion  from  this  example  is  that  for  a  normal 
beam,  with  well  defined  physical  properties,  the  method  is  capable 
of  giving  good  results  and  quite  high  accuracy  can  be  obtained  for 
the  first  ftw  modes  by  using  a  20-lumped-mass  model, 
b.  Hollow  Fectangular  Box-Beam 

In  his  experiments  on  the  whipping  induced  by  underwater 
explosions,  Chertock  (G3)  used  a  thin-wailed,  stiffened,  rectangular 
box-beam.  The  dimensions  are  shown  in  Figure  G1  .  In  this  case  the 
structure  cannot  be  defined  quite  so  exactly  as  for  the  solid  beam. 
The  top  of  the  box  had  a  4"  wide  opening  along  its  entire  length, 
and  during  tests  a  cover  was  bolted  on  to  close  the  opening. 

The  degree  of  fixity  of  the  cover,  and  its  strength,  will  affect  the 
cross-sectional  inertia  considerably.  It  is  assumed  that  the 
bolts  were  really  tight  and  that  the  cover  was  simply  equivalent 
to  the  missing  material.  Ignoring  the  stiffeners,  the  elastic 
data  is 

E  =  29.10^  psi,  v  =  0.3,  4=6  inches 
I  =  9.51  in4,  Ag  =  .529  in2. 

The  total  ballasted  weight  was  105  lb  so  that  the  draft  was  2.70 

g 

in.  Then  H  =  ^  =  1.668  and  the  added  mass  coefficient  C  =  1 .40 

(again  from  Lewis’  paper  (G2)).  The  lengtfy/breadth  ratio  L/B  = 
13.33  so  that  J  =  0.903  and  this  time  the  effect  of  the  3-D 
flow  correction  is  quite  significant. 


In  the  calculation  of  the  added  musses  the  factor  *£»  is 
included  to  ailow  for  the  free  surface  effect.  The  lumped  mass 
data  is:- 

m.  =  5.25  lb;  m  .  =  8.69  lb;  m  .  =-  5.25  ib;  k.  =  1.947  lb/in 

1  W 1  W  X  1 

r.  =  50  lb. in”;  r  .  -  26.0  lb. in”;  r  .  =  15.8  lb.in^; 

1  wi  '  wi  w 

k  .  =5.85  lb. in 
ri 

Here  the  structure  lumped  inertia  is  a  rough  estimate  since  the 
exact  figure  would  vary  from  mass  to  mass  (the  stiffeners  being 
discrete  ana  several  transverse  bulkheads  being  present). 

A  table  of  results  for  the  same  cases  as  for  the  uniform 
beam  is  reacily  constructed: - 

Cases:  1;  'Exact'  uniform  beam  frequencies  (without  shear) 

2;  Calculated  excluding  shear  and  rotairy  inertia 

3;  Calculated  excluding  shear  but  including  rotary  inertia 

4;  Calculated  including  shear  but  excluding  rotary 
inertia 

5;  Calculated  including  both  shear  and  rotary  inertia 
using  the  Taylor  formulation  of  shear  area. 

5a;  Calculated  including  both  shear  and  rotary  inertia 
using  the  Taylor  formulation  of  shear  area. 

6;  Measured  frequencies 
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FREQUENCIES  FOR  CHERTOCK’S  BOX  BEAM  MODEL 


»A 

0 

0 

86 

238 

466 

771 

2A 

0 

0 

88 

244 

479 

794 

3A 

0 

0 

87 

235 

451 

726 

4A 

o 

0 

86 

227 

414 

628 

5A 

0 

0 

85 

221 

399 

602 

5(a)A 

0 

0 

85 

217 

389 

580 

6A 

. 

0 

0 

81 

191-197 

- 

- 

1W 

- 

- 

53 

146 

286 

473 

2W 

1.07 

1  .26 

54 

150 

294 

487 

3W 

0.96 

l  .34 

54 

146 

284 

462 

4 * 

1 .16 

1.17 

53 

139 

254 

385 

5W 

i.ll 

1  .22 

53 

137 

249 

376 

5(a)W 

1.11 

1  .22 

52 

135 

242 

362 

6w 

1.15- 
•  1.35 

U15- 

1.35 

49- 
'  51 

98- 

108 

- 

- 

The  agreement  this  time  is  much  less  satisfactory.  The  ’exact’ 


uniform  beam  results  still  agree  well  with  the  appropriate  case  3» 
and  the  remaining  calculated  results  are  all  quite  consistent. 
Unfortunately  there  is  a  considerable  difference  between  the 
calculated  results  and  those  measured  on  the  model.  For  the 
fundamental  modes  in  air  and  water  the  calculated  frequencies  are 
higher  than  those  measured  by  4.7%  and  4%  respectively,  and  for  the 
three  node  modes  the  corresponding  figures  are  12%  and  31%.  Since 
the  lumped  mass  model  has  been  shown  to  be  accurate  for  such  low 
modes,  and  still  agrees  with  the  ’exact’  uniform  beam  solution, 
the  difference  roust  be  due  either  to  non- beam-1 ike  behaviour  of 
the  model  or  to  incorrect  assignment  to  parameter  values  (masses, 
sectional  inertias,  shear  areaseic.). 
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In  addition  to  the  natural  frequencies,  the  mode  shapes  and 
longitudinal  strain  distributions  were  also  measured  in  Chertock's 
experiments  and  the  strains  deduced  from  the  curvatures  in  the 
normal  modes  were  found  to  be  rather  greater  than  the  measured 
strains  suggesting  that  either  large  shear  deflections  were 
occurring  or  the  sectional  inertias  were  lower  thanexpected. 

The  above  calculated  results,  in  appropriate  cases,  contain  the 
effect  of  shear  but  it  was  thought  that  perhaps  the  simple  beam 
theory  shear  area  formula 

A 

s  m 

(where  I  is  the  cross-sectional  inertia,  t  the  thickness  of 

material  on  the  neutral  axis  and  m  is  the  moment  of  the  area  on 

one  side  of  the  neutral  axis)  might  not  be  very  accurate  for  the 

thin-walled  box  beam.  By  using  an  alternative  formula,  due  to  Taylor  (G4), 

for  the  shear  area  of  thin  walled  girder,  a  smaller  value  for  the 

shear  area  was  obtained  and  this  gave  the  results  for  case 5(a)  above. 

~  common  practice  in  ship  calculations  is  to  assume  the  shear  area 
is  simply  equal  to  the  area  of  vertical  plating.  The  three  types 
of  shear  area  for  the  box  beam  model  are:- 

2 

actual  cross-sectional  area  =  1 .44  in 

2 

(  area  of  vertical  plate  =  0.576  in 

|  2 
shear  (  beam  type  shear  area  =  0.529  in 

areas  (  2 

(  Taylor  shear  area  =  0.429  in 

It  is  clear  from  the  results  however  that  the  differences  in  shear 
area  do  not  account  for  the  discrepancy  in  frequencies.  An 
alternative  partial  explanation  could  bo  that  the  cross-sectional 
inertia,  I  was  too  large.  As  noted  earlier,  the  model  in  fact 
had  a  4  inch  wide  gap  in  the  centre  of  the  top  plating,  which  was 
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closed  off  during  the  tests  by  ;i  bolted-on  lid.  Assuming  that  the 
lid  was  totally  ineffective,  the  moment  of  inertia  and  shear  area 
become 

I  =  7.51  in^»  Ag  =  0.513  in2 

2 

instead  of  1  =  9.51,  Ag  =  O.529  in  as  with  the  lid. 

The  change  in  inertia  will  reduce  the  calculated  frequencies  by  the 
factor 

1. 

'2 

so  that  the  new  frequencies  in  air  and  water  will  be  approximately 

No.  of  nodes  2345 

frequency  m  ^  346  5i5 

air 

frequency  in  ^  ,  20  215  322 

water 

The  fundamental  mode  frequency  is  now  too  low  but  that  of  the  3  node  mode  is  about 
right.  It  is  possible  that  any  relative  motion  between  the  lid  and 
the  main  box  could  affect  the  modes  rather  differently  and  this 
seems  a  plausible  reason  for  the  differences  in  frequencies. 

This  beam  shows  up  the  usual  features  of  ship  structures: 
the  added  water  mass  is  very  large,  1.6  times  the  structure  mass, 
the  three  dimension  fluid  flow  correction  is  appreciable?  1($  of  the 
added  mass,  and  shear  deflections  are  very  important,  particularly 
for  the  higher  modes.  The  shear  correction  for  the  5-node  mode 
frequency  is  26?o.  The  data  for  this  model  however  are  not  adequate 
to  differentiate  between  the  3-D  flow  corrections  of  Taylor  and  Lewis 
or  between  the  normal  beam-type  shear  area  and  Tqylor’s  shear  area. 


c.  Full  Scale  Ship 


The  data  for  a  destroyer  were  recently  calculated  as  part  of 
a  vibration  test  and  are  given  in  the  following  table.  As  for  the 
unifora  beams,  the  ship  was  divided  into  20  equal  length  sections 
and  the  mass  and  rotary  inertia  of  each  section  was  lumped  at  its 
centre.  The  lumped  masses  are  then  connected  by  19  equal  length 
weightless  beams  whose  elastic  characteristics  correspond  to  crass- 
sections  mid-way  between  the  masses.  If  x  =  o  corresponds  to  the 
bow  and  L  is  the  length  of  the  ship,  the  lumped  masses  are  at  the 
positions 

x.  =  (i  -  ^)L/20  ;  i  =  1,  ...»  20 

and  the  cross-sections  characterising  the  weightless  beans  are  at 
xi  =  iL/20  ;  i  =  1,  19 

The  data  require  some  explanation.  The  ship  mass  and  buoyancy 
are  fairly  easily  obtained  from  the  ship’s  design  drawings.  The 
added  mass  was  obtained  from  Lewis'  formula  by  calculating  the 
parameters  a  and  b  of  the  Lewis  transformation  with  the  correct 
bean/ draft  ratio  and  passing  through  a  selected  point  of  the  cross- 
section.  Several  other  points  on  the  cross-section  were  then 
computed  to  verify  the  acceptability  of  the  fit.  The  appropriate 
values  of  C  and  C  were  computed  and  comparison  of  the  buoyancy 
calculated  from  C  with  the  ship  design  value  gave  further  general 
confirmation  of  the  fit. 
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The  method  is  considerably  more  complex  than  applying 
Prohaska's  curves  but  ensures  a  good  representation  of  the  section 

and,  once  programed  for  a  computer,  is  straightforward.  Choosing 
the  collocation  point  either  at  l/3rd  of  the  draft  above  the  keel, 
or  at  the  point  of  attachment  of  the  bilge  keel  gave  the  most  reason¬ 
able  results.  The  existence  of  a  bilge  keel  in  the  midship  region 
is  .itself  a  complication.  To  allow  for  its  effect  references 
C13  and  C5  both  recommend  the  use  of  Wendel's  correction  factor  (D7)» 
Applying  Wendel's  figures  for  the  square  gives  an  increase  of 

(since  the  beam  is  19  feet  and  the  keel  width  2  feet).  Applying 
the  approximation  of  section  D  (the  draft  was  only  12  feet)  gives 
the  value  14$.  Since  both  of  these  values  appeared  to  be  rather 
large  an  alternative  approach  was  tried.  The  radius  of  curvature 
at  the  turn  of  bilge  was  rather  larger  than  the  width  of  the 
bilge  keel  so  it  was  assumed  that  the  flow  near  the  keel  would 
be  fairly  similar  to  the  flow  along  a  flat  wall  with  a  plate  project¬ 
ing  out  at  right  angles.  The  latter  flow  is  simply  one  half  of 
the  two-dimensional  flow  normal  to  an  isolated  flat  plate.  This 
flow  is  well  known.  The  velocity  of  the  flow  past  the  plate  car. 
be  related  t.o  the  vertical  velocity  of  the  ship  section  by  computing 
the  velocity  which  would  exist  at  the  root  of  the  bilge  keel 
according  to  Lewis'  transformation,  it  the  keel  were  not  present. 
Suitable  resolution  of  the  resulting  force  on  the  keel  then  gives 
an  appropriate  added  mass.  The  method  is  described  fully  in 
Appendix  1.  The  correction  for  the  bilge  keel  is  found  by  this 
method  to  be  between  Jjo  and  4$  and  this  is  the  value  which  is 
included  in  the  table.  In  fact  the  keel  is  not  particularly 
important  in  whipping  calculations  since  using  the  full  Wendel 
value  instead  of  the  value  used  here  would  reduce  the  basic  2-node 
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vibration  frequency  by  about  2^  only.  For  vibration  purposes  this 
difference  might  be  more  significant. 

To  correct  the  added  mass  for  the  three-dimensional  flow  effects, 
the  Taylor  correction  factor  for  the  2-node  mode  of  vibration 
was  used.  For  this  ship,  with  a  length/bean  ratio  of  9,3  the  factor 
is  0,807  and  so  is  again  quite  significant;  without  it  the  added 
mass  would  be  750  tons  greater  and  the  basic  frequency  would  be 
dropped  by  around  Lewis*  reduction  factor  for  this  case  is 

0.847  which  would  increase  the  added  mass  by  about  150  tons  and  so 
reduce  the  frequency  by  a  little  over  1?£. 

Calculation  of  the  rotary  inertias  is  extremely  time  consuming 
and  was  not  carried  out  for  this  ship.  It  had  however  been  carried 
out  earlier  for  a  fairly  similar  ship  and  the  values  for  that  ship 
have  been  scaled  to  apply  to  the  present  one. 

The  cross-sectional  moments  of  inertia  for  the  beam  sections  were 
calculated  in  the  usual  way  with  all  continuous  longitudinal  material 
up  to  and  including  the  main  deck  being  considered  effective.  The 
forecastle  deck  was  also  included  entirely  right  back  to  the  break 
of  fo'c'sle.  The  deck  houses  however  were  so  short  that  no 
correction  was  necessary.  One  common  method  of  estimating  the 
shear  area  is  to  take  the  area  of  all  vertical  plating.  This  was 
tried  but  led  to  frequencies  which  were  too  high.  The  shear  areas 
calculated  by  beam  theory,  i.e. 

A  =  — 
s  m 

ere  about  half  the  values  given  by  the  vertical  plating  method 
and  lead  to  much  better  frequencies.  Even  this  method  however  is 
questionable  since  it  depends  heavily  on  b,  tne  thickness  of  the 
side  plating  at  the  neutral  axis.  For  a  few  cross-sections  the 
neutral  axis  occurs  within  inches  of  the  waterline,  where,  for  araour 


reasons,  there  is  an  abrupt  change  in  plating  thickness.  Applying 
the  above  equation  rigorously  could  lead  to  the  shear  area  changing 
by  a  factor  2  if  the  neutral  axis  »oved  a  few  inches.  Clearly  so*e 
kind  of  average  is  required.  Taylor's  method  for  the  shear  areas 
of  thin  walled  box  girders  applies  well  in  this  case  and  several  of 
the  inertias  calculated  by  the  normal  beam  formula  were  recalculated 
in  this  way.  The  results  confirmed  the  beam  formula  exceptionally 
well. 

The  normal  mode  frequencies  given  by  the  computer  program  for 
the  four  cases: - 

1.  Shear  deflections  suppressed,  rotary  inertia  suppressed 

2.  Shear  deflections  suppressed,  rotary  inertia  included 

3.  Shear  deflection  included,  rotary  inertia  suppressed 

4.  Shear  deflection  included,  rotary  inertia  included 

are  shown  in  the  table,  together  with  the  actual  frequency  where  this 
was  known. 

CALCULATED  DESTROYER  FREQUENCIES 


Frequency  (Hz) 

Number  of  Nodes 

Calculated 

Measured 

1 

2 

3 

__4 

0  Heave 

0.21 

0.21 

0.21 

0.21 

- 

1  Pitch 

0.23 

0.23 

0.23 

0.23 

- 

2  ) 

1 .68 

1.66 

1.54 

1.53 

1.59 

3  )c 

3.98 

3.93 

3.13 

3.1 1 

3.20 

)  ° 

'  <0  TJ 

6.90 

6.75 

4.67 

4.64 

4.77 

5  £  § 

>2 

1 1 .02 

IO.64 

6.41 

6.37 

6.29 

6  ) 

7  > 

16.04 

15.59 

8.09 

8.05 

- 

22.54 

21.13 

9.61 

9.56 

- 

The  agreement  is  quite  satisfactory  for  the  four  vibration 
nodes  for  which  an  experimental  frequency  is  known,  except  for  the 
basic  Mde.  Even  this  mode  is  quite  adequate  for  whipping 
calculations  where  only  a  few  cycles  are  of  interest.  Once  again 
rotary  inertia  is  a  very  minor  effect  whereas  shear  is  very  important, 
reducing  even  the  basic  mode  by  9,a  and  for  the  higher,  6  and  7  modes, 
the  frequency  is  more  than  halved.  The  calculated  and  measured 
mode  shapes  are  compared  in  Figure  2.  Tne  agreement  is  again  very 
good  except  for  the  basic  mode.  This  shape  however  could  not  be 
compared  directly  as  the  vibration  generator  was  incapable  of 
providing  a  sinusoidal  input  at  such  a  low  frequency  and  an.estioated 
contamination  of  the  first  mode  by  the  second  has  been  removed 
from  the  plotted  experimental  values.  Part  of  the  frequency 
discrepancy  may  be  due  to  errors  in  the  assumed  mass  distribution. 

A  sister  ship  in  a  similar  displacement  condition  was  found  to  have 
a  frequency  of  1.51  Hz. 
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Explosion  Induced  Forcing  Function 

Following  Chertock’s  original  analysis,  which  was  found 
satisfactory  for  the  models,  the  compressible  phase  of  the 
explosion  pressure  wave,  and  the  migration  of  the  bubble  of 
explosion  products,  are  both  neglected  here.  They  will  be 
discussed  later.  Under  these  assumptions  the  flow,  in  the 
absence  of  the  ship  and  any  boundary  surfaces  (eg  the  free  surface 
or  the  bottom)  will  be  spherically  symmetric  and  if  the  volume  of 
the  bubble  is  V(t)  the  velocity  at  a  distance  r  from  the  bubble 
centre  will  be  given  by 

U  =  i/^2  <*» 

The  acceleration  of  the  water  will  therefore  be 


.  y  2V  dr  V  V 

u  =  - 5 - 5  —  =  - j - j-* 

4*r  4xr  dt  4*r  8x  r 

The  second  term  clearly  decays  much  more  rapidly  than  the  first 

away  from  the  bubble  and, in  any  case,it  is  consistent  with  the 

earlier  development  of  the  forces  via  strip  theory  to  neglect 

2 

forces  of  the  type  ^pu",  so  the  second  term  is  neglected  here. 
The  possible  magnitude  of  this  and  similar  terms  is  investigated 
later.  The  water  acceleration  is  therefore,  in  the  present 
approximation. 


Referring  to  the  geometry  shown  in  Figure  Hi,  this  acceleration, 
when  resolved  vertically  at  the  ship  axis  gives 


a  cost  =  2  ;  R.2  -  D2  +  H2  +  (x.-x  )2 

1  1  A *R-  i  4*R-  1  1  * 

l  l 
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To  allow  for  the  presence  of  a  free  surface  at  the  ship  axis, 
an  image  bubble  of  opposite  sig:.  to  the  real  one  must  be  added  at 
the  reflection  position  above  the  free  surface.  At  the  free 
surface  the  image  will  have  the  effect  of  doubling  the  vertical 
accelerations  and  cancelling  out  any  horizontal  ones,  so  that  the 
final  vertical  acceleration  is 


Substituting  this  into  equation  (F6)  gives  the  modal  equation 


"  2  ,  ;;  i)  k-i  (rawk  +  PAke*yk/Kk 

ai  *  "i  ai  =  xiv ;  -  s  — » - r —  >H4) 

(m  ♦  m  )y. 

k=l  k  wk  ki 

for  the  case  where  rotary  inertias  are  neglected.  The  same 
equation  applies,  out  with  an  obvious  modification  to  when 

rotary  inertias  are  included.  The  similarity  of  equation  (H4) 
to  equation  (BlO),  given  by  Chertock's  analysis,  is  clear.  Tnc 
factor  2  between  the  expressions  is  due  to  including  the  effect 
of  the  free  surface  in  equation  (4)*  Equation  (H4)  is  slightly 
less  accurate  than  (BlO),  due  to  the  ’lumping*  procedure,  but  is 
more  general  since  it  is  not  restricted  to  ’proportional  bodies’. 
It  is  however  still  restricted  to  slender  bodies  where  strip 
theory  is  applicable. 

The  equations  (H4)  retire  only  initial  conditions  and 
specification  of  the  bubble  volimc  history  for  their  solution. 

The  equations  governing  tne  pulsations  of  a  non-migrating  bubuu. 
are  well  known  except  veiy  late  in  the  bubble  collapse  .hose  where 
mathematical  instabilities  m  tie  equations  reflect  an  observed 
physical  instability.  An  excellent  account  of  the  equation-  .mi 


their  derivation  is  given  by  Cole  (ill),  hut  a  brief  outline  is 
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included  here  as  it  provides  a  useful  lead-in  to  the  modifications 
for  migration  effects. 

The  velocity  in  the  water  is  given  try  equation  (Hi)*  i.e. 

•  2* 

V  a  a 


where  a  is  the  bubble  radius.  The  kinetic  energy  is  therefore 


ri»i,  2.2,  .  3.2 

K.E  =  /  ^pu  .4^r.dr  =  2xpa  a 

a 

Tlie  potential  energy  of  the  bubble  is  a  combination  of  the 
internal  energy  of  the  gas  products  and  the  potential  energy  of 
the  hole  in  the  water,  so  that 


T-  1 


-  7ta  poo 


where  p  is  the  gas  pressure,  ani  p  the  ambient  pressure  at  the 

OO 

bubble  centre.  Since  the  pulsations  are  very  rapid  the  adiabatic 
equation  of  state  is  appropriate,  i.e. 


k  Ky 

p  =  = 


where  W  is  the  charge  weight.  If  Eq  is  the  total  energy  tied 
up  in  the  bubble  motion  and  Z is  the  hydrostatic  head  of  water 
at  the  bubble  (£2),  the  equation  of  energy  conservation  is 


E„  ■  “  * PE  •  *  5  ^‘3gZo  *  —  ;  Y-I  3(r-i> 

(Y-1)  <f )  a 

This  cay  be  written  more  simply  in  terms  of  the  non-dimensional 
variables  x  =  a/L,  t  =  t/T; 
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Equation (H6)  is  the  differential  equation  satisfied  by  the  bubble 
radius  and  applies  so  long  as  the  bubble  is  spherical  and  remote 
from  all  surfaces.  In  practice  small  bubbles  fairly  deep  remain 
spherical  throughout  their  motion  except  for  very  short  periods 
centred  about  the  bubble  minimum.  At  these  times  the  bubble  surface 
is  unstable  and  greatly  distorts.  The  bubble  minimum  and  maximum 
radii  are  given  by  (H6)  when  x  =  o,  i.e. 


v3  k  1 

’  ?lY=tT‘ 


and  for  small  k  are  approximately 

i  r  i 


k3T7H)  ,  ,  kY-t/3fo-o  . 


minimum 


^  1-  j  !  1  +  (y-  j)k  ;  maximum 


(H7) 


(H8) 


The  various  constants  in  these  equations  have  been  determined 
empirically  by  Arons&Cj.  He  found  that,  for  TNT,  the  correct 
variation  of  bubble  period  (from  explosion  to  minimum)  with  depth 
could  only  be  reproduced  by  equation  (H6)  if 


03 


Y  ~  5/4 

and  then  the  actual  value  of  the  period  is  correct  if 

E  =  eW 
o 

where  e,  the  energy  per  unit  mass  of  charge  is  approximately  490 
cal/gm.  Comparing  further  the  experimental  peak  bubble  pressure 
with  the  predicted  value  gives 


k  .0552  Z 


Together  with  the  theoretical  expression  for  k  given  after  equation 
(H6)  this  provides  the  gas  equation  of  state 


P  =  8.l5pJ* 

where  p  is  in  kbars  and  p  in  gras/cc.  Jones  (H3)  had  some  time 
earlier  predicted  theoretically  that  this  equation  should  be 


P  =  7.8p'c 


for  the  low  pressures  involved  during  most  of  the  bubble  period. 

With  Aron's  values,  the  scale  constants  L  and  T  and  the 
bubble  maximum  radius  a  and  period  T  are  given  approximately  by 

S3  p 

L  =  15.6  w’^/zj^3  ft 


(HI  0) 


T  =  2.94  W,/3/Z05/6  seconds 


12.6  W1//3/Z  ,//3  ft 


T  r  4.36  W1//3/Z  5^6  seconds 

P  ° 


where  W  is  the  weight  of  the  churge  in  lbs  of  TNT  and  Z  is  the 

o 


hydrostatic  head  in  feet  of  water.  The  bubble  volume  acceleration 
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is  then  given  by 


V(t)  =  ^  ^  (x3)  =  1220  W'^Z  2/3(x3)  ft3/ sec2 

3  T‘  ° 

••5 

where  x  ,  by  differentiating  (6),  is  given  by 


(H11) 


x3  =  3x2  x  „  bxx2  =  ^  l/x3  -  4  -  (3r-4)k/x3^j 


(HI  2) 


The  pressures  in  the  water  around  the  bubble  are  also 

proportional  to  x'ror,  except  close  to  the  bubble  where  -’^pu  terms 

are  important,  the  pressure  is  simply 
* 

P  =  P$ 


Here  (/>  is  the  velocity  potential  for  the  flow  around  the  bubble  i.e. 


<P  = 


V 

4xr 


so  that 


and  with  the  same  units  as  before 


(H13) 


pr  =  £-  V  =  1  .341  W’/V^V)  Psi-ft.  (H14) 

ifT,  O 

The  equations  (III 0)  for  the  period  and  amplitude  of  the  bubble 

motion  apply  strictly  only  to  TNT.  To  produce  bubble  details  for 

close  comparison  with  Chertock's  experiments  the  basic  bubble 

equations  can  be  used.  The  detonator  which  produced  Chertock's 

bubble  was  fired  13  inches  below  the  model,  i.e.  at  a  depth  of 

16  inches  below  the  surface  so  that  Z  =  33  +  1.33  =  34*33. 

o 

The  observed  radius  was  5.4  inches  and  the  period  23.5  msecs. 

This  period  will  have  been  reduced  somewhat  from  the  free-water 
period  given  by  equation  (HI 6)  due  to  the  proximity  of  the  wster 
surface.  The  empirical  correction  factor  for  this  effect  is 
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a 

(see  reference  H4)  (l  -  0.2  -if)  i.e.  0.932  for  this  bubble. 

u 

The  equivalent  free-water  period  is  consequently  23.5/.932  =  25.2 
msecs.  Now 


T  T  x  ,  l/2  -r 

_£  T  _ra  / _ 3\  _ m 

a  =  L  x  "  l2gZ  ’  x 
m  m  0  o  m 


(H15) 


so  that  the  ratio  of  period  to  amplitude  is  independent  of  the  energy 

release  and  depends  only  on  Z  and  on  k  (through  -r  ,  x  )•  From 

0  b  m 


reference  H3,  the  solutions  of  equation  (H6)  for  y  *  1 .25  give 


k 

X 

m 

X 

0 

T 

m 

vx« 

0 

1.0000 

0.0000 

1.492 

1.492 

0.1 

0.9645 

0.0464 

1  .4821 

1.537 

0.2 

0.9235 

0.1172 

1.4804 

1 .603 

0.3 

0.8742 

0.2031 

1  .491 8 

1.706 

0.4 

0.8100 

0.3064 

1 .5120 

1.867 

From  equation  (Hl5),  for  Chertock’s  detonator, 
a  2gZ  l/o  T 

^  =  (^r-V/2  -E  =  1.520 

ra  m 

so  that  k  =0.066  (for  TNT  at  the  same  depth  the  value  would  have 

been  0.134)  and  x  =  0.978,  t  =  1.483.  Then,  to  give  the  correct 
m  m 

radius. 


L  =  5.52  inches  and  so  T  =  1 6.95  msecs. 
With  these  values. 


pr  =  V  =  r?  x3  =  5  pLgZ  x3  =  18.71  x3  psi-inches 
**  J  4*'T  9  0 


These  results  apply  only  up  to  the  first  bubble  minimum  since  near 
the  minimum  energy  is  lost  from  the  pulsation  both  by  acoustic 
radiation  and  by  cooling  from  water  spray  projected  into  the  bubble 


as  its  surface  becomes  unstable.  By  a  similar  process  to  the  above, 

again  making  allowance  for  the  surface  effect,  the  observed  values 

17*6  msecs,  3.9  inches  and  15.2  msecs,  5.2  inches  for  the  period 

and  radius  for  the  second  and  third  pulses  lead  to  the  values: 

k  =  0.120,  T  =  13.1  msecs,  L  =  4.25  inches  (second  pulsation) 

k  ;  0.217,  T  •  to .7  msecs,  L  =  3.48  inches  (third  pulsation) 

These  give  results  for  the  first  three  bubble  pulsation  which 

are  compared  with  the  measured  values  of  pr  in  Figure  2.  Chertock's 

theoretical  values  for  pr  were  very  similar  to  these  but  were 

obtained  for  a  constant  k  =  0.134  (corresponding  to  TNT)  and  neglecting 

the  surface  effect.  His  length  and  time  scales  were  arbitrarily 

chosen  to  give  the  observed  values.  The  steady  increase  in  k  found 

here  for  the  later  pulsations  is  to  be  expected  since,  from  equation 

(H6),  k  is  proportional  to  (E  )~^  and  for  successive  pulses  E 

o  o 

becomes  steadily  snaller. 

To  produce  pressure  or  particle  velocity  histories  similar  to 
those  in  Figure  2  the  differential  equation  (H6)  has  to  be  integrated 
and  then  (Hi 2)  gives  the  necessary  values  of  ?.  The  integration 
is  most  conveniently  carried  out  using  the  technique  described 
in  reference  (H6) .  This  technicpie  uses  Tschebyschef f  polynomials 
to  evaluate  the  definite  integral  which  equation  (H6)  gives  for 
the  time  taken  to  expand  from  minimum  radius  to  maximum  radius. 

The  technique  may  be  readily  extended  to  give  the  indefinite 
integral  for  ”c(x).  The  method  gives  high  accuracy  with  very  few 
terms. 
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Description  of  Computer  Program  I 

To  simplify  the  application  of  the  preceeding  equations  tc  actual 
cases  a  computer  program  was  written  in  FORTRAN  and  implemented  on 
Atlas,  KDF-9  and  IBM  7090  computers.  Copies  of  the  program  are  held 
at  the  Naval  Construction  Research  Establishment  and  a  detailed 
specification  for  the  program  structure,  the  data  format  and  units 
and  types  of  output  which  can  be  produced  are  given  in  Reference 
[il],  A  brief  description  of  the  program  is  given  here.  It  is 
possible  that  the  normal  mode  part  of  the  program  could  have 
application  to  ship  vertical  vibration  problems.  The  program 
produced  all  the  results  described  in  Section  G. 

The  elastic  nature  of  the  ship  structure  is  specified  by  giving 
the  number  n  of  lumped  masses  to  be  used  in  the  representation 
(up  to  4®  on  Atlas),  the  cross-sectional  inertias  and  shear  areas 
of  the  (n-t)  beams  connecting  the  masses  and  the  elastic  constants 
(Young's  Modulus  and  Poisson's  Ratio)  of  the  hull  material.  The 
ship  inertial  properties  are  specified  by  giving  a  ship  mass  to  be 
associated  with  each  lumped  mass  together  with  an  added  water  mass 
and  the  mass  of  water  displaced  by  the  section  of  ship  associated 
with  the  lumped  mass.  If  rotary  inertia  is  thought  to  be 
significant  the  rotary  inertia  to  be  associated  with  each  lumped 
mass, due  to  the  ship  mass  must  also  be  given.  Added  water  rotary 
inertias  are  assigned  by  the  program.  Normal  pitch  and  heave 
motions  are  included  and  at  each  lumped  mass  an  immersion  force 
per  unit  submergence  may  be  given. 

From  the  elastic  structure  data  the  computer  builds  up  the 
stiffness  matrices  in  equation  (C4)  and,  if  rotary  inertia  is  not 
required,  eliminates  the  rotations  as  in  equation  (C6).  The 
inertia  mass  data  are  then  used  tc  build  up  the  mass  matrices  in 
equation  (E9).  The  transforation  described  in  section  F  is 


used  to  produce  a  symmetric  matrix  whose  eigenvalues  and  vectors 
are  found  by  anHGW  type  of  eigenvalue  subroutine. 

If  a  series  of  positions  (up  to  20  on  Atlas)  are  specified 
along  the  ship  the  program  then  calculates  the  coefficients  y^(x), 

IL(x)  and  S^(x)  from  equation  (F8)  for  use  in  equations  (F7) • 

The  program  then  integrates  the  modal  coefficient  equation  (H4)  and 
recombines  the  equations  as  in  (F7)  to  produce  histories  of  the 
relevant  quantities  at  each  of  the  specified  positions  aloi.g  the 
ship.  It  is  necessary  to  give  the  position  of  the  charge  relative 
to  the  ship  in  the  data  so  that  the  program  can  form  the  forcing 
function  coefficients  in  equations  (H4).  The  bubble  volume 
history  is  specified  by  a  subroutine.  A  table  of  values  of 

(t,  x3(t))  is  used  for  interpolation. 

Normally  the  program  carries  out  the  integration  assuming  that 
initially  the  ou  and  their  time  derivatives  ^  arc  all  zero  but, 

if  a  particular  initial  velocity  distribution  is  specified  in  the 
data,  the  program  n solves  it  into  an  equivalent  set  of  non-zero 
a^'s  to  start  the  integration. 

In  order  to  integrate  the  differential  equations  by  step-by- 
step  methods  a  time  increment  is  required.  For  the  linear 
oscillator  eqrations  the  Runge-Kutta  integration  routine  is  stable 
if  this  increment  is  less  than  about  l/6th  of  the  period  of  the 
highest  frequency  mode  present.  The  data  needed  for  the  program 
includes  the  number  of  steps  to  be  used  in  each  such  period. 

Ten  steps  seems  to  give  quite  accurate  results  judged  by  comparing 

with  a  similar  integration  using  20  steps  per  period.  Although 

this  increment  ensures  stability  of  the  differential  equations 
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it  is  far  too  large  to  give  an  adequate  representation  of  the 
bubble  forcing  function  near  the  positive  pressure  'spikes*.  In 
these  regions  the  program  automatically  reduces  the  step  length  to 
ensure  the  representation  is  reasonable. 
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Complete  Example 

For  the  box  model  described  in  section  G,  together  with  the 
bubble  forcing  function  derived  in  section  H,  the  computer  program 
was  used  to  predict  the  motion  of  the  model  in  Chertock's  experimental 
tests.  Figure  Jt  compares  the  measured  and  predicted  result  for 
the  two-node  mode  coefficient  for  a  charge  fired  at  l6  inches  beneath 
the  free  surface,  directly  under  the  centre  of  the  model.  The 
agreement  for  the  early  motion  is  quite  good  but  in  the  later 
motion  the  experimental  amplitudes  become  slightly  less  than  predicted 
and  also  the  two  curves  gradually  swing  out  of  phase.  The 
differences  in  phase  are  possibly  due  to  a  reduction  in  the  model 
frequency  due  to  the  amplitude  of  the  motion.  The  difference 
first  shows  significantly  at  around  27  mseconds,  at  which  time  the 
strains  in  the  model  are  about  3/4  of  yield.  Prior  to  this  time, 
the  difference  between  the  two  curves  is  very  small  and  explainable 
by  the  2Hz  spread  in  the  measured  frequency.  The  reduced  experimental 
amplitudes  could  well  be  due  either  to  (Lmping,  which  is  neglected 
in  the  theoretical  model,  or  to  minor  inaccuracies  in  the  bubble 
forcing  function,  which  is  not  known  at  all  accurately  in  the  later 
stages. 

FigureJl  also  shows  the  lowest  frequency  rigid  body  mode  and 
it  is  clear  that  the  vertical  displacement  of  the  model  largely 
follows  the  changes  in  the  bubble  volume.  The  mode  shape  shown 
is  not  actually  either  s  pure  heaving  or  pure  pitching  mode.  The 
frequencies  of  the  two  motions  were  so  close  that,  to  the  accuracy 
specified  in  the  eigenvalue  routine  in  the  computer  program,  they 
were  identical  and,  although  the  routine  always  produces  orthogonal 
vectors  for  multiple  eigenvalues,  the  first  is  arbitrarily  chosen 
in  the  corresponding  eigenvector  space.  For  this  geomelry  the 
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two  rigid  body  modes,  when  taken  together,  produce  a  pure  heaving 
motion. 

Figures  J2,  J3  and  J4  show  the  relative  effect  of  variations 

in  the  charge  position:  longitudinally,  athwartship  and  in  depth, 

for  the  first  three  whipping  modes.  For  the  2-node  mode  the 

experimental  results  given  by  Chertock  have  been  added  to  the 

f^g<_res  and  the  agreement  is  excellent.  The  effect  of  mode  shape 

on  the  longitudinal  variation  is  evident.  For  the  transverse  and 

depth  variations  no  curves  have  been  given  for  the  second  (3-node) 

whipping  mode  since  at  all  points  in  the  transverse  plane  through 

the  model  centre  (where  these  results  apply)  the  amplitude  induced 

for  this  mode  is  zero,  the  centre  being  a  node.  From  the  curves 

it  may  be  seen  that  the  amplitude  of  the  third  mode  decays  slightly 

more  rapidly  than  that  for  the  basic  mode,  both  transversely  and 

with  depth.  To  produce  such  comparisons  it  is  not  necessary  to 

carry  out  complete  integrations.  Since  the  equations  are  linear 

and,  except  for  very  minor  variations  in  the  bubble  period  with 

depth,  the  bubble  volume  acceleration  is  the  same  for  all  charge 

location^  the  modal  amplitudes  are  proportional  to  the  forcing 

function  coefficients  (X.  in  section  H).  It  is  these  coefficients 

1 

which  are  compared  in  Figures  J2,  J3  and  J4» 

The  results  of  this  example  sb,--  that  the  lumped  mass/weightless 
beam  idealisation  of  the  ship,  together  with  the  strip  theory 
approximation  for  fluid  effects,  can  represent  the  whipping  motions 
of  a  small  scale  model  quite  accurately.  Comparison  with  Chertock’s 
original  work  shows  that  the  present  theoretical  model  is  as 
accurate  as  his  original  formulation,  while  being  slightly  more 
general  in  that  it  applies  also  to  non-proportional  bodies.  It. 
is  also  more  adaptable  and  can  readily  be  extended  to  include  other 
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effects  such  as  gravity  migration  and  plastic  deformation  of  the 
target.  The  use  here  of  strip  theory  shows  up  particularly  clearly 
the  nature  of  the  forces  acting  on  the  ship  and  their  distribution 
alorio  U . 


RELATIVE  AMPLITUDE  RELATIVE  AMPLITUDE  RELATIVE  AMPLITUDE 


i  rw yzi 


-0-6 


BASIC  MODE 


-  THEORETICAL 

o  0  o  EXPERIMENTAL 


DISTANCE  FROM  CENTRELINE  -INCHES 


DISTANCE  FROM  CENTRELINE  -  INCHES 


SECOND  MODE 


\  DISTANCE  FROM  CENTRELINE  -  INCHES 
40  \  50  80  IOO  120 


THIRD  MODE 


RELATIVE  VARIATION  OF  WH-°PING  MODE  AMPLITUDES 


WITH  LONGITUDINAL  CH#  RGE  POSITION 


FIGURE  J  2 


79 


iTIVE  AMPLITUDE  RELATIVE  AMPLITUDE 

_ A. 


TRANSVERSE  DISTANCE  FROM  CENTRELINE  -  INCHES 


TRANSVERSE  DISTANCE  FROM  CENTRELINE  -  INCHES 


RELATIVE  VARIATION  OF  WHIPPING  MODE  AMPLITUDES 
WITH  TRANSVERSE  CHARGE  POSITION 


FIGURE  J  3 


79 


16  24  32 

»  ■  * 

DEPTH  BELOW  FREE  SURFACE  -  INCHES 


THIRD  MODE 


In  the  development  of  the  theory  as  presented  in  sections  E, 

F  and  H,  it  has  been  assumed  that  the  motion  of  the  ship  is  independent 
of  the  initial  shock  wave,  that  in  fact  the  ship  behaves  as  though 
initially  a  highly  compressed  bubble  of  gas  were  suddenly  released 
from  rest  to  accelerate  and  expand  outward  and  then  pulsate*  In  fact 
the  bubble  at  no  stage  can  be  considered  to  be  at  rest  and  the 
water  surrounding  it  does  not  accelerate  to  high  velocities  in  the 
way  this  model  suggests.  What  actually  happens  is  that  when  the 
detonation  wave  in  the  explosive  reaches  the  water  it  transmits  a 
shock  wave  into  the  water  and  almost  instantaneously  raised  the 
bubble  surface  to  a  very  high  velocity.  Similarly  in  the  surrounding 
water  the  shock  wave  raises  the  velocity  to  a  high  value  and  all 
these  initial  high  velocities  gradully  decay  until  they  fall  to  the 
incompressible  velocities,  after  which  normal  incompressible  bubble 
motion  takes  place.  Even  the  initial  bubble  size  assumed  in 
section  II  is  fictitious,  being  u.-ually  considerably  larger  than  the 
original  charge  size.  It  is  merely  the  size  of  bubble  which,  if 
the  water  were  really  incompressible,  would  pulsate  with  the 
period  and  maximum  radius  which  are  observed  in  experiments. 

Chertock's  original  work  (31)  gave  a  qualitative  argument, 

based  on  an  application' of  Green’s  theorem  to  the  flow  potentials, 

^  • 

assuming  them  to  be  solutions  of  the  wave  equation,  why  it  should 
be  expected  that  the  acoustic  snock  wave  phase  would  not  produce 
significant  motion.  In  fact  his  arguments  ao  not  rigorously 
apply  to  the  type  of  discontinuous  function  involved  in  the  early 
interaction  and,  in  addition,  did  not  really  include  all  the  inter¬ 
action  waves  (reflected,  diffracted,  etc.).  In  addition,  for 
surface  ships  at  ieast,  non-linear  effects  are  certainly  present 
during  the  initial  phases  as  cavitation  effects  are  well  established 


(Kl)  during  the  shock  wave  hull  plating  interaction  phase  of  the 
motion.  Cavitation  has  the  effect  of  largely  preventing  aqy 
reflection  of  incident  momentum,  so  ensuring  that  the  momentum  in  the 
shock  wave  is  fully  absorbed  by  the  ship. 

In  consequence  the  neglect  of  the  shock  wave,  or  at  least  its 
replacement  by  a  zero  velocity,  high  pressure  bubble,  must  be 
considered  to  be  justified  principally  by  the  experimental  evidence. 
On  small  scale  the  agreement  between  the  observed  and  predicted 
model  motions  certainly  fully  justifies  the  assumed  artificial 
bubble.  However,  the  scaling  laws  of  underwater  explosions  are 
such  that  it  is  generally  true  that  shock  wave  effects  are  relatively 
much  smaller  than  bubble  effects  on  small  scale.  The  assumption 
cannot  therefore  be  considered  established  for  a  larger  scale 
by  the  small  scale  experiment. 

There  are  however  other  heuristic  reasons  for  expecting  that 
it  should  still  be  valid  at  least  for  some  conditions.  For  a 
target  deep  enough  for  cavitation  not  to  be  significant  Murray, 

(K2),  showed  by  an  exact  solution  that  when  a  spherical  acoustic  wave 
interacts  with  an  infinite  neutrally  buoyant,  circular  cylinder 
(considered  to  be  rigid  in  cross-section  but  infinitely  flexible 
in  shear  along  its  length),  the  transverse  velocity  of  the  cylinder 
asymptotically  tends  to  the  nearby  water  particle  velocity  of  the 
incident  wave.  He  found  in  fact  that  initially  the  cylinder 
jumped  to  a  considerably  higher  velocity  which  then  decayed  to 
the  incompressible  value  after  a  period  equivalent  to  several 
transit  times  of  the  cylinder  diameter  by  the  acoustic  wave. 

For  non-neutral ly  buoyant  circular  cylinders,  the  associated 
two-dimensional  problem  of  a  plane  wave  incident  on  a  circular 
cylinder  of  elastic  cross-section  has  been  considered  by  Haywood 
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(K3)  who  found  that  in  this  case  the  cylinder,  after  again 
reaching  a  fairly  high  initial  velocity,  slowed  down  to  the 
asymptotic  velocity 

- (Kl ) 

1  +  p 

where  u  is  the  undisturbed  water  particle  velocity  at  the  cylinder 
centre,  due  to  the  passage  of  tne  shockwave.  In  this  equation  a 
is  the  cylinder  radius  and  p  is  the  mass  of  the  cylinder  per  unit 
lengtn.  The  previous  incompressible  assumption  gives  the  same 
result,  as  the  fluid  impulse  which  produces  a  velocity  u  in  the 
water  imparts  a  momentum  (p^  +  u^)u  Per  unit  length  of  the 

cylinder.  This  momentum  is  given  to  the  cylinder  and  its  added 
water  mass  so  that,  if  the  resulting  cylinder  velocity  is  v, 

v(|i  +  pw)  =  (pw  +  Pw)u.  (K2) 

For  a  circular  cylinder  the  added  mass  per  unit  length,  pw,  is 
7tpa“  as  also  is  the  buoyancy  per  unit  length,  p^.  Whence  0t2)  gives 
the  same  value  for  v  as(Kl). 

These  results  suggest  that,  at  least  when  cavitation  is  not 
important,  the  incompressible  assumption  at  the  start  will  lead  to 
quite  good  results  in  the  later  motion.  Where  cavitation  leads 
to  most  of  the  momentum  in  the  shock  wave  being  absorbed  by  the 
target  again  the  deviation  from  the  incompressible  theory  should 
be  small  because  the  total  momentum  in  the  initial  positive  phase 
is  generally  at  least  8QJo  due  to  the  ’incompressible'  pressures 
which  follow  behind  the  very  short  duration  compressible  shock 
wave  head  of  the  positive  phase.  The  subsequent  motiom  of 
the  ship  and  its  entrained  water  will  therefore  again  be  largely 
as  though  the  incompressible  motion  only  were  present  at  the 


start 
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An  alternative  to  the  Strip  Flow  Method  for  the  Added  Water  Mass 


The  usual  theory  for  calculating  the  hydrodynamic  forces  on  a 
ship,  applied  as  described  in  Section  D,  is  clearly  an  extremely  simple, 
and,  for  the  2-node  mode  at  least,  very  accurate  method  of  determining 
such  forces.  For  the  higher  whipping  modes  the  experimental  evidence 
is  very  limited.  The  strip  theory  does  however  have  cne  possible 
drawback  in  the  present  application.  In  the  development  of  the 
normal  mode  form  of  the  equations  described  in  section  F,  all 
the  »ode  frequencies  and  shapes  were  found  simultaneously  through  the 
eigen  values  of  a  single  matrix.  This  matrix  naturally  involves 
the  strip  theory  added  water  masses  and,  through  them,  the  correction 
factor  J  which  compensates  for  three  dimension  flow  effects.  In 
consequence,  in  the  development  used  here  the  same  value  for  J  is 
applied  to  all  the  modes.  This  conflicts  with  the  customary  use 
of  the  factor  ,  where  a  different  value  for  J  is  used  for  each  aide 
shape  (including  the  pitch  and  heave  modes).  To  use  the  customary 
method  here  would  mean  developing  a  separate  matrix  (using  the 
appropriate  value  for  J)  for  each  mode,  and  extracting  just  the 
appropriate  eigen  root  and  vector  from  each  matrix.  In  the  examples 
given  earlier  the  simpler  alternative  was  adopted  of  applying  the 
2-node  mode  correction  factor  to  the  added  masses  (since  that  is 
the  mode  contributing  most  to  flexural  stresses)  and  ignoring 
the  fact  ttet  all  the  other  modes,  pitch  and  heave  included, 
will  be  somewhat  in  error. 

The  type  of  correction  factor  which  should  be  used  in  vibration 
calculations  has  been  the  subject  of  a  great  deal  of  discussion  ever 
since  strip  theory  was  first  used.  The  two  independent  oiiginators 
of  the  theory,  Lewis  (Li)  and  Tailor  (L2)  both  presented  reduction 
factors  but,  due  to  assuming  slightly  different  boundary  conditions 
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they  arrived  at  different  reduction  factors.  Which  factor  to  use, 
Lewi-.'*  or  Taylor's,  isstill  discussed  today.  A  furtt.:*-  difficulty, 
which  has  been  raised  in  numerous  papers  ((L3)  contains  a  useful 
discussion)  is  that  it  has  never  been  rigorously  established  that 
the  correction  factor  should  be  constant  over  the  length  of  the  ship, 
as  is  implicitly  assumed  in  normal  useage.  In  fact,  in  at  least 
one  recent  paper  (L4)  a  correction  factor  has  been  presented  which 
does  vary  along  the  length,  although  in  this  particular  case  no 
basis  for  the  given  variation  was  presented. 

Because  of  these  basic  difficulties  in  the  customary  strip  theory, 

particularly  as  applied  to  the  whipping  problem,  the  alternative 

approach  described  here  lias  been  developed.  The  fundamental  weakness 

of  strip  theory  is  that  the  force  it  predicts  on  each  cross-section 

is  related  only  to  the  motion  of  that  cross-section  and  not  to  the 

probably  slightly  different  motions  of  neighbouring  cross-sections. 

In  the  notations  used  in  section  E,  the  vertical  hydrodynamic  force 

S  .  on  the  ith  mass  is,  according  to  strip  theory,  S  =  -  a  y. 

W1  wi  wi i 

and  the  total  vector  of  hydrodynamic  forces  is 

K  *  ■  o-i) 

where  is  a  diagonal  added  mass  matrix.  This  standard  fora  of 
strip  theory  <quation  clearly  cannot  be  valid  for  general  motions. 

If  one  section  of  the  ship  is  accelerating  upward  it  will  certainly 
experience  a  downward  hydrodynamic  force  so  a  diagonal  term  is 
certainly  expected.  However,  its  upward  acceleration  will  induce 
a  downward  fluid  acceleration  around  it  as  fluid  moves  to  fill  the 
space  vacated  beneath  it.  Part  of  this  flow  will  be  around 
neighbouring  cross-sections  and  so  will  induce  a  downward  forces 
on  these  section^  too.  For  the  above  equation  to  be  generally 
true  then  it  is  necessary  for  to  contain  also  off-diagonal 


elements.  In  fact  all  elementi-  of  Mw  should  be  nonzero  but,  since 

each  section  will  affect  principally  its  nearest  neighbours,  the 
magnitude  of  the  off  diagonal  elements  will  decay  rapidly  away  from 
the  main  diagonal.  The  problem  is  to  approximate  the  eier\nts  of 
in  a  simple,  readily  calculable  manner.  Sine-  the  general  ship 

problem  is  much  more  complex,  attention  is  directed  first  to  the  case 
of  an  exisymmetric  ship.  A  method  for  applying  the  axi-syametdc 
results  to  a  normal  ship  is  described  later. 

Mathematical  Model  tor  an  Axi-Syametric  Ship 

For  such  a  ship  it  is  possible  to  satisfy  fairly  well  the  fluid 
.-oundary  condition  of  equality  of  hull  and  fluid  velocities  along 
the  normals  to  the  hull  by  means  of  a  vertical  dipole  distribution 
along  the  axis  of  the  ship.  This  distribution  should  normally  be 
continuous  but,  if  the  ship  is  considered  to  be  divided  into  a 
number  of  sections  and  the  number  is  large  enough,  then  it  should 
be  reasonable  to  assume  that  the  line  distribution  in  each  section, 
i.  has  a  constant  strength 

Since  the  velocity  potential  at  the  point  (r,S,x)  due  to  a 
dipole  at  the  point  (o,o,s)  (referred  to  the  cylindrical  polar 
coordinates  of  Figure  1)  is  given  by 


<t>  = 


ir  cos  6  ,  0 . 

vy$  »  whcre  u  = 


(r  +(x-s)  ]- 

the  potential  due  to  a  line  distribution  of  dipoles  of  strength 
jju  per  unit  length  and  extending  from  Xf  \/2 


is 


<^i(r,Q,xl=  fr.rcosQ 


f  XW2 

j 

Xi-l/2 
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(L2) 


whence 


The  radial  velocity  u  produced  by  the  distribution  is  given  by 

r 


"xi1**  ■  -  3T  *  i*ireose  f 


l  [r  +(x-x. 


i+l/2' 


r  2  ,  X2-.3/2  J 

[r +(x-xi_1/2)  J  J 


The  total  velocities  produced  by  all  the  distributions  are  therefore 


n  n 

u  (x)  =  2  u  .(x)  and  u  (x)  =  Z  u  .(x) 
r  j=J  rj  x  j=1  xj 

The  boundary  condition  at  the  hull  is  that  the  velocities  of 
the  fluid  and  the  hull,  resolved  in  the  direction  of  the  normal  to 
the  hull,  should  be  equal*  For  the  axi-syautetric  shape  considered, 
the  angle  X  between  the  normal  to  the  hull  and  the  radial  direction 
is  given  by  X  =  -  tan  *b#  where  b(x)  is  the  radius  of  the  ship,  and 


b'=  Ub/dx.  The  boundary  condition  is  therefore, 

vcosX  -  v'bsinX  =  u  cosX  +  u  sinX  (L5) 

r  x 

where  v(x)  is  the  distribution  of  vertical  velocity  along  the  length. 
This  condition  allows  for  both  shearing  and  flexing  of  thediip. 

For  the  more  interesting,  lower  frequency,  modes  of  vibration, 
v'b  will  be  small  compared  to  v  since  the  wavelength  will  be  much 
greater  than  the  half  beam.  X  will  also  normally  be  small  except 
possibly  in  the  immediate  neighbourhood  of  the  stern.  The  term 
v'bsinX  will  therefore  be  very  small.  ux»  the  flow  along  the 
ship,  is  produced  partly  by  the  variation  in  v  along  the  ship  (i.e. 
by  v')  and  partly  by  the  changing  shape  of  the  ship  (b').  It  too 
will  therefore,  normally  be  small  when  v'and  b'  are  small.  It  is 
generally  possible  then  to  simplify  the  boundary  condition (15)  to 

v  =  ur  (L6) 

except  when  rapid  changes  occur  in  either  the  ship  underwater  cross- 
sectional  area  or  in  the  velocity  distribution  along  its  length 
(e.g.  for  high  modes).  This  simplification  in  the  boundary 
condition  reduces  slightly  the  amount  of  data  needed  to  specify 
the  ship  and  its  motion.  It  is  roughly  equivalent  to  neglecting 
rotary  inertia  in  the  dynamic  equations  of  the  ship  itself. 

Neither  of  the  boundary  conditions  (L5)  or  (L6)  can  be 
satisfied  everywhere  by  the  assumed  velocity  potential.  However, 
either  can  be  satisfied  at  up  to  n  'collocation'  points  along  the 
ship  and  if  these  points  are  suitably  chosen  the  boundary  condition 
will  not  be  seriously  violated  anywhere.  The  most  convenient 
choice  for  the  collocation  points  is  at  the  mid-points,  x^,  of 
the  ship  sections  and  satisfying  equation(L6)  at  these  points  gives 
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(y  ~  u)  =  ~  A(i 


(L7) 


where  b  is  the  maximum  value  of  b(x),  y(x)  the  vertical  velocity  of 
the  ship,  u(x)  the  vertical  velocity  at  the  ship  due  to  the 
explosion  and  the  matrix  A  =  (a .  .)  where 

*  J 


.  -  f  x.-x.+e/2  X.-X  .-&/2 

L_  i  J _  _  _ JL-  J _ 

'Aj  b^  1  [b.%(x.-x.  +e/2  )2]1/2  [b2i+(x >*] 


5TT72  + 


bi2(xi~x  aI/2)  b  i(xi~x  ,-l/2)  •» 

[b  Mx  -x  fl/2)  2]5/2  [b2.  +  (x.-x  V/2  )  2]^2  J 

X  1  J  x  *  J 


U8) 


Equation (L7)  may  then  be  inverted  to  determine  the  strengths 
H  of  the  line  distributions  necessary  to  satisfy  the  boundary 

condition  (L6)  for  the  given  velocity  distribution  y  .  This  gives 


p  =  b2AH  (y  -  u) 


The  upward  vertical  force  per  unit  length  on  the  section  at 
the  point  x^^  is  given  by 


f.  =  -2/  p  cos  0  b.d© 
1  o  1 


where  p  =  p<£  is  the  fluid  pressure, 

n 


Here  the  potential  d>  =  d>  +  2  <#> .  since  the  fluid  motion  due  to 

0  j=] 

the  explosion  (and  having  the  potential  </>q)  is  included.  Thus 

f.  =  -2pJ*  h<*>  cos  0d0  =  2Pa.b.2/*  cos20d0-2a  2n  a  A.  /1tcos20d0 
1  r  o  iT  rix'o  T.  l  J=1  ij  j  o 

=  Pi/j  <L’> 
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x.-x-.-e/  2 

1  J  1 

] 

In  the  vicinity  of  the  ship,  <f>Q  has  been  approximated  by  <pQ  «  -u^  = 
-ub^  cos  6  on  the  ship  surface.  The  velocity  u  is  the  explosion 
induced  velocity  at  the  ship  axis. 

Since  in  any  actual  motion  conditions  along  the  ship  will  be 
continuous,  the  upward  forcc/unit  length  at  the  centre  of  the  section 
will  approximate  the  average  force/unit  length  over  the  whole  section 
and  the  total  upward  force  S^,  on  the  section  will  be  approximately 


ib.\(x  -x.+t/2  n1'2 

i  l  j 


>.  =  &  f .  =  xpb.24  -  xpb  £n  p.  .u. 

l  l  r  i  i  r  .  .  rijri 

J  “  * 


"  •  •  •  t  * 

i.e.  ^  =  Mwu  -  *pbBp  =  *wu  ~  *plrBA”  (£  -  u)  (LIO) 

*v  ^ 

where  B  is  the  matrix  (j3.  .)  and  M  is  the  diagonal  matrix  with 

1 J  * 

elements  mwi=xpb^  l  .  The  matrix  is  exactly  the  same  ’buoyancy' 
matrix  as  found  before  in  section  E. 

The  retired  inertial  water  mass  matrix  is  therefore  given  by 

Mw=*b3pBAH  (LI  1) 

This  matrix  depends  only  on  the  shape  of  the  ship  and  is  completely 
independent  of  the  type  of  motion  (or  vibration  mode).  With 
known  the  force  distribution  for  <*ny  vertical  acceleration  distribution 
y  is  readily  found.  A  short  computer  routine  has  been  written  to 

compute  the  non-dimensional  matrix  BA~*  and  the  force  distributions 


deduced  from  it,  by  equation  (LtO),  for  given  distributions  of 
vertical  acceleration. 

Coup ari son  of  Results  with  known  Exact  Force  Distributions 

There  are  two  foras  of  fluid  flow,  suitable  for  comparison,  lor 
which  exact  solutions  are  known:  -  the  flows  around  a  vibrating 
prolate  spheroid  and  around  an  infinite  circular  cylinder  whose 
transverse  velocity  varies  sinusoidally  along  its  length.  The 
exact  results  for  both  cases  have  been  compared  with  the  results 
from  the  foregoing  analysis, 
a.  Vibrating  Prolate  Spheroid 

The  methods  involved  in  the  solution  of  the  flow  about  a 
prolate  spheroid  are  discussed  in  some  detail  by  Lamb  (L5). 

Lewis  (Li)  gave  the  first  solution  in  connection  with  transverse 
shear  vibrations  and  Taylor  (L2)  gave  a  second  solution  using  a 
different,  more  realistic  boundary  condition  involving  flexure 
as  well  as  shear.  A  more  general  solution  by  Macagno  and 
Landweber  (L6)  contains  both  the  earlier  solutions  as  special  cases. 

The  coordinate  system  used  in  these  solutions  is  the 
ellipsoidal  set  (£,p,*>)  illustrated  in  Figure  Li.  The  curves  of 
constant  £  and  p  represent  confocal  systems  of  ellipsoids  and 
hyperboloids  of  two  sheets  respectively.  In  terms  of  the 
cylindrical  coordinates  (r,0,x) 

r  =  k(C2-1),/2  (1 W/2 


x  =  k£p 

If  £=  CQ  represents  the  given  ellipsoid,  assumed  to  be  of  length 
2a  and  diameter  2b,  then  k  and  £q  are  given  by 


k  =  (a2-bV/2  =  a(l-pV/2 


where  £  =  t/a 


«..<«  -pV,/2 

If  the  vertical  velocity  distribution  along  the  axis  of  the  ellipsoid 
is  v(»/a) 

then  the  boundary  condition  on  the  surface  of  the  ellipsoid  will  be 

u  =  -  7^  =  (v  cos  X  -  ~  b3in  X)  cos  9 

n  dnrr  dx 

^ o 

which  is  the  type  of  condition  used  by  Taylor,  allowing  for  flexure 
as  well  as  shear.  If  v(x/a)  is  a  polynomial  of  degree  N,  then  the 
solution  for  the  potential  can  be  written 


<f>  =  £  a  d  (C)lJ  (p)  cos 


where  (p),  (O  are  associated  Legendre  functions  of  the 

first  and  second  kinds  respectively.  Substituting  this  expression 
for  £  into  the  boundary’  condition  gives 


jjvM  -  SM  ajr  :  <(*•=' 

J=C0 

(since  x  =  ap  on  the  ellipsoid  C=C0)»  By  assumption  v(j/a) 
is  a  polynomial  so  that 


v  =  2  v  (i) 

«  n  a 
n=l 


and  the  boundary  equation  becomes 


N+l  dP  dQ  N+1  o  „  , 

2  “du  k  dT  *  *n^  =  2  )P  3vnH  » 

n=1  ^  K  ^  C=C  "  n=1  n 


which  determines  the  coefficients  aR  of  the  potential  in  terms  of  the 


coefficients  v  of  the  velocity  distribution.  This  equation  is 
most  easily  solved  by  first  integrating  with  respect  to  p,  giving 


N+1  N+1 

n=1  n=1 


-a  § 

dQn 
n  I 

ank 

d£  -i 

to  +1 

gives, 

3  xn 

where  represents  -a^jj-  *  Multiplying  by  P^p)  and 


for  Legendre  Polynomials 


N+ 1  *  /  .  \ q2 

X  =  —  v  1  f  (m=\9..9$R+\) 

w  n  n  mn 


where  lmn  =  /  p  PB(p)dp  =  o,  n  <  m 


or  (n-«)  odd 


2  n(n~l )  •  •  •  (n-m+2) 
*rn+ro+1  )(n+m-1 )  ...(n-on-3) 


(n-m)  even 


From  these  equations,  given  v  ,  x  is  easily  found  and  so  then. 

n  m 


from  its  definition,  is  a  • 

m 

With  the  coefficients  a^  known  in  the  expression  for  the 

velocity  potential,  the  force  distribution  on  the  ellipsoid  is 
readily  found  from 


f (x)  =  -2  J  b(x)p  cos  0d6  where  p  =  ptf> 


i.e. 


N+1 


f(x)  =  -xpb(x)  ^  anQi*tt0>Pl'(Va) 


(LI  2) 


of 


A  .snort  computer  routine  has  been  written  to  compute  the  vaiues 
f(x)Apb2  given  b/a  and  values  for  the  coefficients 


It  is  limited  to  velocity  distributions  which  can  be  represented  by 
quartic  equations  but  this  is  sufficient  for  the  lowest  few  modes. 

Figures  L2  and  L3  compare  the  exact  force  distributions  given 
by  the  above  equation  with  the  results  given  by  the  approximate 
3-D  flow  method,  for  two  ellipsoids  with  b/a  ratios  of  l/5th  and 
l/lOth.  The  values  given  by  normal  strip  theory,  i.e. 


N+1 


=  (l-x2/a2)  2  v  (i)n_1 


n  a 


(LI  3) 


xpb  n=l 

have  also  been  added  to  the  figures.  The  Lewis  correction  factors 
have  been  used  for  the  strip  theory  results  since  these  are  known 
for  four  of  the  five  nodes.  The  values  used  for  the  coefficients 
Vj»«.*vj  for  these  cases  were 


>v»«>^co  e  f  f  ic  i  ent 

Mode 

• 

V1 

• 

V2 

• 

V3 

• 

V4 

• 

v5 

Heave 

1 

0 

0 

0 

0 

Fitch 

0 

1 

0 

0 

0 

2-node  vertical 

-0.2 

0 

1 

0 

0 

3-node  vertical 

0 

-0.429 

0 

1 

0 

4-node  vertical 

0.0274 

0 

-0.534 

0 

1 

the  first  four  of  these  distributions  being  those  used  by  Lewis. 

The  fifth  distribution  has  nodes  at  0.1  55L,  0.38L,  0.62L,  0.845L, 

L  being  the  total  length.  These  nodes  are  in  about  the  correct 
position  for  a  4-node  ship  vibration  mode  but  the  relative  amplitudes 
are  not  good  since  the  distribution  gives  too  much  prominence  to 
the  ends  and  too  little  to  the  central  section. 
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for  the  i/a  =  l/lOth  ellipsoid  (for  a  ship  with  a  length  to  beam 
ratio  of  10)  the  results  of  the  3-D  approximation  and  strip  theory 
noth  agree  well  with  the  exact  analysis  and  there  is  little  to 
choose  between  them.  Near  the  centre  of  the  ellipsoid,  where 
changes  in  the  radius  are  smallest,  the  3-D  approximation  is  slightly 
better  but  at  the  ends,  where  the  radius  is  changing  rapidly,  strip 
theory  is  the  better.  For  the  b/a  =  l/5th  (L/B  =  5)  ellipsoid  the 
results  are  similar  but  the  divergence  of  both  approximations  from 
the  exact  solution  is  quite  serious  near  the  ends,  even  for  modes 
as  low  as  the  second  vibration  mode.  Again  there  is  little  to 
choose  between  the  approximate  3-D  flow  and  strip  theory. 

However,  in  each  case  the  strip  theory  uses  a  different  correction 
factor  for  each  mode  whereas  the  3-D  approximation  deduces  all  the 
forces  from  a  single  matrix  equation. 

One  point  which  should  be  borne  in  mind  in  these  comparisons 
is  that  the  boundary  conditions  are  not  quite  the  same  for  the 
three  methods.  Strip  theory  necessarily  uses  a  shear  type 
boundary  condition  since  it  can  only  relate  the  force  on  a  section 
to  the  vertical  velocity  there.  The  3-D  approximate  analysis 
too,  based  on  equation  (L9)»  also  includes  only  shearing  motions. 

The  definition  of  matrix  A  (equation  L8)  is  based  on  the  approximate 
form  of  the  boundary  condition  equation  (L6),  which  is  equivalent 
to  including  shear  only.  The  exact  solution  includes  flexure 
as  well  as  shear.  The  necessary  modifications  to  the  3-D  approxi¬ 
mation  in  order  to  satisfy  the  more  exact  boundary  equation  (L5)» 
and  the  results  it  gives,  are  described  later. 


b.  Infinite  Circular  Cylinder 


Taylor  (L2)  was  the  first  to  consider  this  example  in  a  ship 
vibration  context.  The  infinitely  long  cylinder  is  assumed  to  have 

a  transverse  velocity  distribution 

v (x)  =  v  cos  kx 
o 

The  velocity  potential  for  this  distribution  Taylor  found  to  be 

v  K  (kr) 

<p  _  JLJ -  cos  6  cos  kx 

kfc'(kb) 

where  b  is  the  cylinder  radius  and  is  a  modified  Bessel  function 

of  the  second  kind.  The  force  distribution  m£*y  be  found  as  before 
and  is 

_.  2  vK  (kb) 

f(x)  =  -2  /*p  cos  0bd6  =  -2pb/  <f> cos0d0  =  -  xpb  -2 -  coskx 

o  o  kbK'(kb) 

The  wavelength  of  this  distribution  is  of  course  X  =  2 r/k. 

The  example  is  not  ideal  for  comparing  with  the  3-D  approximation 
since  the  latter  is  firmly  tied  to  a  body  of  finite  length  (although 
it  could  obviously  be  modified  to  deal  with  a  periodic  distribution 
of  infinite  length).  To  compare  results  the  approximate  analysis 
was  applied  to  a  finite  length  of  cylinder,  3  full  wavelengths  being 
the  selected  length  in  the  hope  that  the  central  wavelength  would 
be  free  of  end  effects.  Since  the  computer  program  for  the  3-D 
approximation  is  restricted  to  only  20  sections  this  gives  only 
6  different  dipole  distribution  regions  per  wavelength.  The 
results  are  compared  with  the  exact  values  from  the  above  equation 
in  the  table.  In  this  case,  because  the  cylinder  is  of  uniform 
diameter  there  is  no  difference  in  the  two  forms  of  boundary 
condition. 


(LI  4) 
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\*A 

r/b 

V2 

5*/b 

7Vb 

3*/b 

1  lx/b 

13Vb 

5V2 

17Vb 

S  exact 

mm 

0 

.448 

.448 

0 

-.448 

-.448 

0 

.448 

3-D 

WSSSk 

-.00 

.445 

-.0014 

-.0014 

-.449 

-.452 

-.0140 

.383 

-.631 

0 

.631 

.631 

0 

-.631 

-.631 

0 

wm 

■BaM 

-.609 

-.0003 

.609 

.609 

-.0013 

-.611 

-.616 

-.0189 

Km 

-.762 

0 

.762 

.762 

0 

-.762 

ES9 

0 

.762 

HgQH 

-.741 

-.0001 

.741 

.741 

-.0004 

-.744 

-.0096 

.674 

- - - n 

20  exact 

-.826 

0 

.826 

.826 

0 

-.826 

-.826 

0 

.826 

3-D 

-.839 

-.0007 

.839 

.839 

.0007 

1 

-.639 

-.840 

-0015 

.823 

- - ‘-•y—  -  .  . . . 

Values  or  F(x)/xpb  vq  for  an  Infinite  Cylinder 


For  all  values  of  X/b,  the  forces  are  clearly  very  accurate  near 
the  centre  being  within  of  the  exact  value  in  all  cases.  The  points 
given  were  the  only  ones  used,  indicating  the  coarseness  of  the 
representation;  equivalent  in  a  ship  representation  to  using  twenty 
points  to  represent  the  fifth  vibration  node  (6  nodes).  Since  the 
values  from  the  3-D  approximation  are  constant  to  at  least  x/\  = 

7Vb,  the  differences  from  the  exact  result  are  due  to  the  coarse 
mesh  rather  than  the  finite  length.  For  the  ~h/b  =  5  case  the 
lengths  of  the  sections  are  nearly  equal  to  their  radii,  whereas  for 
the  X/b  =  20  case  the  sections  are  each  nearly  9  radii  long. 

Application  of  the  Technique  to  non-axi symmetric  Ships 

Clearly  it  would  be  possible  in  principle  to  extend  the  3-D 
approximation  by  adding  distributions  of  higher  multipoles  along 
the  ship  axis,  and  determining  their  strength  by  satisfying  a 
boundary  condition  around  the  circumference  of  the  ship  as  well 
as  along  it.  Alternatively  (and  more  accurately)  surface 
distributions  of  sources  could  be  used  over  the  surface  of  the 
ship.  Either  technique  would,  however,  require  very  large  amounts 


i 


of  data  to  represent  the  ship  shape  and  would  require  the  inversion 
of  a  very  large  matrix.  The  success  of  the  strip  method  using  Lewis 
sections  for  the  cross-sectional  shapes  points  to  simpler  approaches. 

In  the  solution  of  the  two  dimensional  flows  about  ship-type 
cross-sections,  although  the  velocity  potential  may  in  fact  consist 
of  a  super-position  of  2-D  multipoles  of  all  orders,  the  added  mass 
of  each  section  depends  only  on  the  dipole  term,  the  shape  of  the 
section  determining  its  strength.  In  the  3-D  case  it  should,  there¬ 
fore,  be  approximately  correct  to  account  for  the  shape  of  the  section 
via  the  strengths  of  the  dipole  distributions. 

Since  Lewis'  original  work,  it  has  been  customary  to  represent 
the  added  mass  per  unit  length  of  the  2-0  cross-sections  in  the  form 


where  b  is  the  half  beam  of  the  section  and  C  is  a  constant  depending 
on  the  section  shape.  Values  for  C  have  been  computed  for  a  great 
variety  of  shapes  (References  LI,  L2,  L7,  L8,  L9,  LlO  and  Ll  1 ) .  The 
above  value  for  the  added  mass  is  also  t'»at  due  to  a  circular 
cylinder  of  radius 


b 


equiv 


=  b  C 


(L15) 


For  a  non-axisymmetric  ship  then,  each  cross-section  may  be  compared 
with  the  known  shapes,  C  determined,  and  then  (15)  gives  the 
appropriate  radius  for  an  axi symmetric  approximation  to  the  actual 
shape.  There  are  no  reasonably  simple  non-axisymaetric  3-D  flows 
with  exact  solutions  which  can  be  used  for  comparison  but  the 
procedure  should  give  reasonable  results.  Certainly  in  the  2-node 
vibration  mode,  in  the  central  section  of  the  ship  where  the  added 
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mass  is  most  important,  the  technique  will  give  very  good  answers 
since  in  this  region  the  flow  is  very  nearly  2-D  and  in  the  2-D  case 
the  method  is  exact* 

Use  of  the  Full  Boundary  Condition 

The  results  given  so  far  for  the  3-D  approximation  have  all 
been  based  on  the  approximate  equation  (L6)  for  the  boundary  equation 
at  the  ship  surface  instead  of  the  correct  equation  (L5)«  If  the 
full  equation  is  to  be  used,  then  the  bending  rotations  y^  of  the 

ship  must  be  included*  If  the  expressions  (L3)  and  (L4)  for  ur(x) 

and  ux(x)  are  substituted  into  the  boundary  condition  (L5)  the 

result  is 


h  -  uoi  ♦  bibi'(fl  -  a'oi>  *  A  Vi 


where 


“ij  *  “ij  *  b2kibI«bi2  *  <*i  -  *}  ‘  V2>*J 
[bi2  *  (xi  -  Xj  -  V2)23*Vi J 


21-3/2 


This  mqy  be  written  in  the  matrix  form 


(j  -  u)  +  B  (r  -  «')  =  ^2  % 

b 


(L16) 


where  B  is  a  diagonal  matrix  with  elements  b.b, u  is  again 

11  ”• 


the  water  motion  at  the  ship  axis  due  to  the  explosion  (calculated 
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-,<rvv^  ,,“,!>-t” .-- - 


as  though  the  ship  were  not  present),  and  u'  is  its  derivative  with 
respect  to  x. 

Inverting  the  equation  gives 


p  =  b2A  '  (y-u)  +  b2A  ’gtf-  u#) 


(LI7) 


and  the  force  distribution  f(x)  follows  as  before.  Since  changes 
in  the  cross-section  are  now  included  it  is  appropriate  in  this 
case  to  deduce  from  f(x)  both  the  force  and  moment  at  each  lumped 
mass  by  means  of  the  integrals 


Xi +1/2  xi+  1/2 

S  =  /f(x)dx-.  M.  =  /  (x-x. )  f  (x)dx 

Xi  -1/2  Xi  -1/2 

•  — 

whence  S  =  M  u  -  xpCp 

<N<  *V 

(L18) 

M  =  R^u'-  xpDjl 


where  the  elements  of  the  matrices  C  and  D  are  given  by 


l.y'l  =  (Hb/2)^  (ru  -r^  *r,_) 


a,+Vr.  a9+^r 

*  nJ?  108  ‘v^  ‘  108  ‘  V^) 


x  .-x  ,+e/2 

where  x^  =  ^  ; 


*1  = 


Vl^V 

l+b.'* 


x.-x  .-e/2 

x„  =  -±-4- 


t 


a.  = 


b.b.'/e  + 

1  1  '  1 


1+b. 

1 


»2 


**1+  =  C(a1+^)2+^]^  ;  r^_  =  [(a^  J$)2+p2]  -  ,  with 
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latifiiyikiiii 


similar  expressions  for  a^,  r2  +  ar,<*  r2-‘ 


diy'«  =  ^2  O  +  b.'  2)  ‘  [ (a.j  +!2)ri+  -  (^->2) r^_  -  (cy-  -£)r2+  +  (a2“^)r2J 


a  +-2+r 


<V-k~ 


♦  - '-TT;  l0«  (a-k.r!j,) - 2 - TiS10*  (a2-k,r~*} 

2(l-b".  *) -1  "l  <,r|-  2(l ib*  V  a2^  2- 


V2al 


V2a2 


+  - <r-n  (r  -  r  ) - 3-7,  (ro  "ro  ) 

(Ub-.  ¥  ,+  1  (i^2)"2  2*  2- 


2  a  2  2  a  2 

c«.  +  ?.  a-  2+^o 

1  x  -  2a  atf2+r  2  x0-2a0  a0+ij+r0 

— JL_i  log  (  i . ,-lt)  +  - 22  log  (2  2±.) 

2\-i  ai  /t.K#2V2  ao^2+ro 


(l+bi  ) 


1 


(Ub.  )■ 


2- 


(L19) 


Thus 


-*»  «■ 
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L  j 

L 

-i 

-1 

CA 

C X 

B 

-1 

-1 

DA 

DA 

B 
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which  gives  the  hydrodynamic  force  and  moment  at  each  lumped  mass. 

Use  of  this  equation  instead  of  equation  E(t)  in  the  derivation  of  the 
equations  of  motion  gives,  in  place  of  equation  E(8),  the  equation 
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Here  not  all  rotary  inertia  terms  have  been  neglected.  The  terms 
affecting  the  displacement  equation  (the  M^2  terms)  are  retained 
since  these  are  the  terms  which  differentiate  between  the  ’shear’ 
and  ’shear  plus  flexure'  boundary  conditions.  Even  when  normal 
rotary  inertia  is  suppressed  these  terras  can  have  some  effect  on  the 
hydrodynamic  forces. 

The  force  distributions  given  by  equation  (Li 8)  (with  u  =  o) 

have  been  added  to  Figures  L2  and  L3  where  they  differ  appreciably 
from  the  simpler  3-D  approach  used  earlier.  The  significant 
improvement  in  the  accuracy  of  the  results  for  the  higher  modes, 
particularly  for  the  smaller  length/beara  ratio  ellipsoid,  is  evident. 
The  results  using  equation  (LI 8)  (based  on  the  full  boundary 
condition  (L5))  are  everywhere  almost  identical  to  the  exact  values. 

In  computing  the  results  from  equation  (LI 8),  in  order  to  ensure 
that  no  more  data  should  be  used  than  is  normally  necessary  in  ship 


vibration  calculations,  thr  slopes  b*of  tne  ellipsoids  were  estimated 
numerically  from  tne  given  radius  tr  at  the  collocation  points.  The 
method  appears  to  be  adequate  for  the  rather  extreme  slopes  involved 
in  ellipsoids  and  so  should  be  quite  satisfactory  for  ship  forms 
where  shape  changes  are  less  severe. 


Modification  of  the  Normal  Mode  Equations 


Because  the  matrix  M  in  the  new  equations  of  motion 


M^y  +  Sy  =  H.,«S  +  O'  (L23) 

»V  *V>  *V 

is  now  no  longer  <1  uia^/nal  matrix,  the  transformation 

z  =  M  y  (L24) 

used  in  section  F  can  no  longer  be  applied.  In  fact  is  no  longer 
a  symmetric  matrix  and  the  transformation  technique  should  not  really 
be  applied  at  all.  An  alternative,  once  asymmetiy  is  accepted  would 
be  to  re-write  the  equation  as 

y  +  M^Sy  =  4 

and  solve  the  homogenous  part 
y  +  'sy  =0 

directly.  This  unfortunately  would  considerably  restrict  the 
available  computer  routines  since  most  of  the  really  efficient  routines 
deal  only  with  symmetric  matrices.  However,  in  ail  cases  so  far 
considered,  the  matrix  M,  has  been  found  to  be  almost  symmetric, 
except  for  elements  corresponding  to  the  end  lumped  masses,  and  there  the 
actual  values  of  the  elements  tend  to  be  small.  To  avoid  having 
to  use  general  matrix  methods  the  matrix  Mj  was  arbitrarily  made 
symmetric  by  replacing  all  the  elements  m^^  by  (m^  +  Wjj^)/2. 


The  forces  resulting  from  this  amended  matrix  were  compared  with 
those  produced  by  the  correct  matrix  for  the  two  ellipsoid  cases 
and  the  two  sets  of  forces  were  found  to  differ  by  no  more  than 
anywhere#  Over  most  of  the  length  of  the  ellipsoid  the  forces 
were  identical.  Since  the  vibration  frequencies  depend  approximately, 
on  the  square  root  of  these  forces,  even  this  small  difference  will 
be  further  reduced  in  the  effect  oi.  the  eigenvalues  and  use  of  the 
modified  matrix  seems  justified. 

Kith  M  symmetric,  a  standard  Choleski  decomposition  routine 
(L14)  may  be  used  to  generate  a  matrix  L  such  that 

LLT  =  (M  +  Mw)  (L25) 

and  the  transformation (L24)  is  replaced  by 


giving 

~  -1  *  -1  * 
z  +  E  z  =  L  M2u  +  L  M,u' 

-1  -T 

where  =  L  SL  ,  a  symmetric  matrix.  The  eigenvalues  and  vectors 

of  Ej  can  be  found  just  as  before  and  the  normal  mode  technique 

can  then  be  applied  exactly  as  in  section  F. 

Application  to  a  Ship  Vibration  Calculation 

For  the  destroyer  considered  in  section  G  the  same  basic  data 

were  used,  with  a  modified  version  of  the  computer  program,  to  predict 

the  normal  mode  vibration  frequencies.  The  results  obtained  are 

compared  in  the  table  with  the  straightforward  strip  theory 

technique.  The  calculated  results  are  for  the  cases 

1.  Normal  strip  theory  with  the  Taylor  correction  for  the  2-node 
mode 


1°5 

tflttiifagf  b  >  i  ■  i  i  i  f 


2,  3-D  approximation  based  on  equation  (Ll  1 )  (the 
approximate  boundary  condition) 

3.  3-D  approximation  based  on  equation  (L22)  (the  full 
boundary  condition). 

In  each  of  these  cases  the  effect  of  shear  deflection  is  included 
but  that  of  rotary  inertia  has  been  suppressed. 


Frequency  (H  ) 
z 

Number  cf 
Nodes 

Calculated 

Heasureo 

1 

2 

3 

0 

Heave 

0.21 

■ 

0.20 

1 

Pitch 

0.23 

0.22 

2 

) 

) 

1.54 

1.51 

1.52 

1.59 

(1.50) 

3 

\  v 

>  1 

)  § 

3.13 

3.14 

3.17 

3.20 

(3.09) 

4 

)  S 

)  is 

4.67 

4.77 

4.81 

4.77 

5 

b.4l 

6.72 

6.78 

6.29 

6 

)  > 

8.09 

8.68 

8.79 

7 

) 

9.61 

10.45 

10.61 

- 

It  is  clear  that  the  method  changes  the  results  little,  which 
is  to  be  expected  since  the  length  to  bean  ratio  of  the  destroyer 
is  close  to  that  of  the  more  slender  of  the  two  ellipsoids,  where 
the  exact,  strip  theory  and  3-D  approximation  force  distributions 
were  found  to  be  all  very  similar.  The  large  difference  between 
the  calculated  results  and  the  raer.sured  values  reflects  mainly  the 
inaccuracy  of  the  ship  data  since  measurements  on  a  sister  ship 
in  a  similar  displacement  condition,  gave  the  values  shown  in 
brackets.  Unfortunately  the  higher  mode  frequencies  were  not 
measured  for  the  second  ship. 

For  a  different  class  of  destroyer  similar  results  were 
obtained  and  are  shown  in  the  following  table, 
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Frequency  (H  ) 
z 

Number  of 
Nodes 

Calculated 

Measured 

> 

2 

3 

0 

Heave 

0.20 

0.20 

- 

1 

Hitch 

0.22 

1  ES  1 

0.23 

- 

2  ) 

1  .64 

1 .64 

1 .66 

1.68 

3! 

Vibration 

3.16 

3.20 

3.22 

3.35 

4  ) 

Modes 

4.92 

4.96 

5.00 

4.97 

5  } 

6.90 

6.96 

7.03 

6.63 

In  each  of  these  cases  the  calculated  mode  shapes  were  virtually 
independent  of  the  hydrodynamic  approximation  used. 

Conclusions  Regarding  the  Hydrodynamic  Flow 

The  results  for  the  force  distribution  on  ellipsoids  show  that 
in  general  the  3-D  flow  approximation,  using  the  full  boundary 
condition,  should  provide  more  accurate  results  than  normal  strip 
theory.  However,  the  difference  is  only  marked  for  small  values 
(=5)  of  the  ratio  (length/beam) .  For  a  value  of  this  ratio  of 
10  the  3-D  flow  approximation  was  still  very  accurate  but  strip 
theory  also  gave  good  results  even  close  to  the  ends. 

The  considerable  difference  between  the  inertial  mass  matrices 
given  by  strip  theory  and  by  the  3-D  approximation,  and  the  close 
agreement  of  the  results  they  give,  suggests  that,  so  far  as 
vibration  of  the  ship  as  a  beam  of  rigid  cross-section  is  concerned, 
little  improvement  from  more  exact  ideal  fluid  theories  is  possible. 

Any  remaining  discrepancies  between  measured  and  predicted  (.ydrodynamic 
forces  are  probably  largely  due  to  non-ideal  flow,  free  surface 
effects,  non-linear  Bernoulli  forces  or  the  distortion  of 
cross-sections. 
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FORCE  PER  UNIT  LENGTH 
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The  theory  as  so  far  presented  has  been  shown  to  give  very 

good  agreement  with  the  limited  experimental  results  available 

but  involves  a  number  of  assumptions,  some  of  them  quite  drastic. 

In  this  section  some  of  these  assumptions  are  discussed  in  more 

detail  to  siio.v  more  clearly  the  limitations  of  the  theory, 

:tost  of  the  approximations  are  those  normally  encountered  in 

snip  vibrations  and  in  consequence  they  have  been  investigated 

extensively  already.  The  limitations  of  beam  theory  for  the 

elastic  representation  of  the  ship  are  beginning  to  be  appreciated 

and  alternative  approaches,  e.g.  (Ml)  have  been  proposed.  In 

general  though,  particularly  for  slender  ships,  beam  theory 

has  been  found  quite  adequate.  In  the  hydrodynamics  of  ship 

2 

vibrations  it  is  customary  to  neglect  the  ’Bernoulli1,  t^pii 
pressure  terms,  and,  for  the  free  surface  condition  even  the 
linear  gravity  terras  arc  neglected.  Both  these  points  have  been 
considered  by  several  authors,  e.g,  (M2)  and  (M3)  and  have  been 
found  small  for  the  vibration  modes  although  they  can  be  significant 
in  the  pitch  and  heave  motions.  For  example,  the  inclusion  of  a 


linear  free  surface  effect  changes  the  added  mass,  and  adds  a 

2 

damping  force  to  the  motion  if  the  'Froudc  number*,  —  (where 

S 

b  is  the  half  beam,  w  the  circular  frequency  and  g  the  acceleration 

of  sravity),  is  less  than  about  4.  For  the  heave  and  pitch  motions 

bw2 

of  the  destroyers  described  earlier,  —  -  1  and  these  modes  and 
the  deduced  rigid  body  motions  may  be  expected  to  be  slightly  in 


error,  but  for  even  the  lowest  frequency  vibration  mode  —  =25 

S 

and  the  approximations  may  be  expected  to  be  very  good. 


Compressible  effects  ir.  ship  vibration  have  been  considered 
by  both  Taylor  (D2)  and  Kaplan  (DS)  and  have  been  found  negligible 


Viscous  effects  seen  generally  to  have  received  iess  attention 
but,  since  with  the  small  amplitudes  involved  in  the  vibration 
frequencies  a  full  wake  cannot  form  beneath  the  shin  when  it  is 
moving  up,  only  skin  drag  is  likely  to  natter.  For  flexural 
vibration  tiie  Reynolds  number  is  of  order  10'  and  such  forces  will 
be  small.  Both  compressibility  and  viscous  effects  lead  principally 
to  damping  which  is  very  small  in  vibrating  ships.  In  ship 
vibration  work  damping  is  usually  attributed  largely  to  internal 
effects  and  structural  damping  rather  than  to  the  hydrodynamics. 

Even  minor  damping  is  important  in  such  work  since  it  iimits  the 
amplitudes  of  vibrations  produced  by  sources  operating  close  to 
one  of  the  natural  frequencies.  In  the  explosion  induced  whipping 
application  however  the  vibration  exciter  (the  explosion  bubble) 
is  only  effective  over  about  two  basic  whipping  cycles  and  the 
effect  of  damping  is  very  small. 

Another  possible  j curce  of  error  in  vibration  work  is  sympathetic 
vibration  of  local  str.iture.  Even  if  the  frequency  of  the  bottom 
of  the  ship  (distorting  as  a  panel,  or  two  panels,  over  its  whole 
width  is  much  higher  than  the  ship  frequency,  it  may  be  forced  in 
and  out  at  the  ship  frequency  and  so  reduce  the  amount  of  water 
being  moved  vdth  the  ship.  The  effect  has  been  investigated  in 
(i./r)  and  shown  to  be  important  for  modes  higher  than  the  4-node, 
mode  but  below  this  the  effect  seems  to  be  small. 

Aside  from  the  general  and  well  known  limitations,  there 
are  several  approximations  which  are  related  to  the  interaction 
between  the  explosion  and  the  ship  and  which  are  consequently 
exclusive  to  the  whipping  problem.  The  principal  such  approximation 
is  tne  determination  of  the  force  (distribution  acting  on  the  ship 
through  the  use  of  the  fluid  accelerations  which  would  be  induced 
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jy  tne  exclusion  at  t.u:  si.  ip  axis  f  tin;  ship  were  not  present, 
in  both  tne  strip  theory  and  the  3-D  flow  tiieory  this  approximation 
implicitly  assumes  that,  in  the  absence  of  the  ship,  the 
acceleration  would  be  uniform  throughout  a  region  around  the  ship 
position  which  is  large  compared  to  the  ship  cross-section.  t'nder 
such  conditions  very'  little  flow  will  be  induced  along  the  ship  axis 
jy  t..e  snip  and  it  can  be  expected  that  the  approximation  will  be  good 
but,  for  fairly  close  charges,  the  spherical  flow  divergence  around 
tr.c  charge  will  certainly  lead  to  non-unifonn  conditions.  The 
forces  will  be  acting  on  a  short  length  of  ship,  considerable  flow 
along  the  axis  will  occur  and  the  undisturbed  flow  at  any  cross-section 
will  certainly  not  be  uniform.  This  soherical  divergence  will 
impose  an  inner  limit  to  the  range  at  which  the  'distant  flow' 
approximation  can  be  expected  to  apply. 

A  second  possible  limitation  of  the  'strip  theory'  type  forces 
is  tuat,  even  for  distant  charges  where  conditions  are  uniform, 
near  the  bow  and  stern  of  the  ship  flow  along  the  axis  (around  the 
ends)  can  be  expected.  The  3-D  approximation  should,  in  principle, 
be  capable  of  allowing  for  such  flow  but  strip  theory  is  not. 

A  third  possible  source  of  error  is  the  neglect  of  Bernoulli 
pressures  uue  to  the  bubble  flow.  Although  such  pressures  can  be 
neglected  when  they  relate  to  the  induced  ship  vibration,  the 
direct  flow  due  to  the  bubble  expansion,  when  the  explosion  is 
close,  could  possibly  lead  to  more  effective  pressures  of  this  type. 

These  three  possible  sources  of  error  are  considered  in  the 
fol iowi ng  paragraphs . 
a.  Spherical  Flow  Divergence 

Because  of  the  linearity  of  the  assumed  forces  on  the  ship, 
the  equations  of  motion  derived  in  section  C  are  equivalent  to  a 
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separation  of  the  problem  into  two  parts:- 

(i)  the  determination  of  the  forces  acting  on  the  ship, 
assumed  rigid  and  motionless,  due  to  the  explosion  bubble. 

(ii)  The  response  of  the  ship  to  these  forces,  calculated  as 
though  the  bubble  were  not  present. 

In  consequence,  to  investigate  the  effect  of  spherical  flow  divergence, 
etc.,  it  is  only  necessary  to  consider  the  forces  on  a  rigid  target  . 

The  most  convenient  idealisation  for  such  a  target  is  an  infinite 
rigid  circular  cylinder.  As  before  the  effect  of  the  explosion  is 
assumed  to  be  equivalent  to  a  point  source  of  varying  strength  e(t), 
c  being  given  by 


e(t) 


V(t) 

43t 


woe re  V  is  the  bubble  volume.  The  problem  to  be  solved  is  then  that 
of  a  point  source  near  an  infinite  circular  cylinder,  the  geometry 
and  coordinate  system  being  as  illustrated  in  Figure  Ml. 

For  the  geometry  shown,  the  source  induced  velocity  at  the 
cynnuer  axis  is 


u  = 


T,1 

(b  +x  ) 


desolving  this  velocity  perpendicularly  to  the  axis  gives 


Since  the  strip  theory  added  mass  per  unit  length  and  ’buoyancy' 

2 

ji  per  unit  length  are  both  simply  Ttpa  for  the  cylinder,  the 
distribution  of  force  per  unit  length,  on  the  basis  of  the  distant 
flow  approximation,  is  given  by 

2. 

Mx>  =  =  2  fb372  =  ~~2~2~572  0*»> 

1  w  w  p  (bWy/z  o2*eV/2 

in  terms  of  the  dimensionless  quantities  g  =  x/a,  =  b/a. 


The  problem  has  been  solved,  exactly  and  independently,  by  at 


least  two  authors,  Savic  (M5)  solved  for  the  potential  flow  in 
order  to  estimate  the  effect  of  the  cylinder  on  the  bubble 
pulsation,  and  Murray  (M6)  solved  the  similar  but  more  difficult 
case  of  the  interaction  between  a  spherical  acoustic  source  and  an 
infinite  cylinder.  The  methods  used  were  also  similar. 

Savic* s  solution  for  the  velocity  potential  is 


Y>  =  2  6n  cos  n  0  / 

n-o  \=c 


oo  K  ( bX ) 

6  /  J! -  [K '(aX)l (rX)-I  '(aX)K  (rX)] 

X=o  K  '(aX)  n  n  n  n 
n 


cos(xX)dX;  r<b 


oo  oo  K  (r?.) 

=  —  2  e  cos  n 0  /  -E -  [K'(aX)l  (bX)  -  I  '(aX)K  (bX)] 

u  n  \  n  n  n  n 

n-o  X=o  KR(aX) 


cos(xX)dX;  r>b 


where  6q  =  1 ,  =  2  (nil)  and  Ir(x),  Kr(x)  are  Bessel  functions  of 

the  second  kind.  On  the  surface  of  the  cylinder,  r  =  a,  the 
Wronskian  relation 


In(x)K'(x)  -  Kn(x)I'n(x)  =  —)/ x 
simplifies  the  potential  to 


(/,  =  -—  2  c  cos  n  8  1  (B,£) 
^  xa  X=o  n  n M 


co  K  (pX) 

where  3 (£,£)  =  /  — - cos  (£X)d X  . 


X=o  XKn'(X) 


Since  the  pressure  is  given  by 


P  =  p<(>. 


the  force  per  unit  length  of  the  cylinder  cun  be  fount!  by 
integrating  round  the  circumference  of  the  cylinder  and  is  given  by 


(pado)coso  =  2ap  <p  cos  0dO  =  2ap  [- 


a 


so  that  the  exact  inertial  force  corresponding  to  (Ml) is 
<x>  K  (pX) 

f.,(x)  =  -Ape  /  -i -  cos  (iX)dX  (li4) 

X=o  XK*(X) 

Figure  M2  shows  the  function  (pX)/\K'(X)  for  several  values  of  p. 

The  curves  show  t’aat  most  of  the  contribution  to  the  integral  will  come 
from  very  small  values  of  X  since  the  function  decays  very  rapidly 
as  X  increases,  except  for  p  very  close  to  unity  (the  charge  very 
ciose  to  the  target).  In  fact 


so  that  the  integral  converges  for  all  (p,£)  except  (3  =  1,  g  =  0, 
where  it  has  a  logarithmic  singularity. 

For  p  >  5  Figure  M2  shows  that  there  will  be  very  little  contribution 
to  the  integral  in  (M4)  from  values  of  X  >  1 ,  For  values  of  X  in  the 
range  0  <  X  <  1,  the  denominator  of  the  integrand,  XKj(X),  is 


represented  quite  well  by  the  first  term  of  its  power  expansion  in  X, 
Thus,  for  p  >  5  the  denominator  may  be  replaced  by  the  leading  terms 
of  its  power  expansion  and  the  integral  can  be  represented  approximately 
by 


f  him 
l  *;m 


cos(gX)dX  —  - 


CJ 

J  XK1  (pX)  cos  (gX)dX  - 
o 


\  J  X  log  X  Kj ( pX)  cos  (gX)  dX  - 


-  J  X3K,  (pX)  cos  (gX)  d\ 

o 


»«,  ,.^2  (V  .  M  -  l.«  2  >3Hi 

-ap  *’«■  2  apa52 


K'h  ere 


00 

Hj  =  J  Kq  (pX)  cos  (gX)  dX  =  V2a  ;  a  =  (p2  +  g2)^ 

o 

H2  =  J  Xlog  X  Kq  (PX)  sin  (gX)  dX  =  -  -S-  hog  ^ j  +  y  -  2  . 
o  2ff  L 

H{  is  a  well  known  integral  given  in  reference  [M7j.  H2  was  deduced 
from  Hj  by  Fraenkel  [M8]  for  the  asymptotic  evaluation  of  a  very 
similar  integral.  Substitution  of  these  egressions  gives  the 


asymptotic  expression,  valid  for  large  p. 


f2(x)  ~  2*f*  , ft2  Pp2,3/2  "  *Pe 

VP  +  t  J 


<f 


P(P‘  +  O'  - 


lOxpe  —  to*  - ■  j  ■»  0  (o  log  o) 

CP2  +  o1/z 

The  leading  term  in  this  expansion  is  simply  the  'distant  flow' 
approximation  given  in  equation  [Ml], 


! 


I 


The  results  given  for  f(x)/4pe  by  the  ’ distant  flow*  approxima¬ 
tion  (equation  [Ml]),  the  asymptotic  expansion  (equation  [M5])  and 
direct  numerical  integration  of  the  exact  integral  in  equation  [M43 
are  compared  in  Figure  M3  and  in  the  following  table. 

Exact  and  Asymptotic  Force  Distributions 

f(x)/4pe 


3.0 

10.0 

!  i  1.972 

!  1.111  i  1.272 
i  1  3.347 

0.776 

0.965 

0.673 

0.390 

0.523 

-0.189 

0.068 

0.053 

0.087 

0.0032 

0.0017 

0.0033 

Exact  (numerical)  (4) 
Distant  flow  (l) 
Asymptotic  (5) 

!  1.243 
l .250  i  i .005 
;  j  2.t 87 

0.710 

0.805 

0.761 

0.381 

0.479 

0.003 

0.069 

0.057 

0.081 

0.0033 

0.0019 

0.0033 

Exact  (4) 

Distant  flow  (l) 

Asymptotic  (5) 

:  |  0.600 

1 .667  |  0.566 
j  0.851 

0.488 

0.497 

0.560 

... 

OOO 

0.070 

0.065 

0.071 

0.0025 

0.0036 

Exact  (4) 

Distant  flow  (l) 

Asymptot ic  (5) 

i  0.246 
2.500  o.25l 
j  0.283 

0.232 

0.237 

0.252 

0.195 

0.201 

0.190 

0.066 

0.066 

0.063 

0.0043 

0.0036 

0.0044 

Exact  (4) 
Distant  flow  (l ) 
Asymptotic  (5) 

,  1  0.060 
'  5.000 1  0.063 

1  0.062 

i _ i _ 

0.060 

0.062 

O.O61 

0.058 

0.059 

0.058 

0.0396 

0.0396 

0.0386 

0.0059 

0.0056 

0.0060 

Exact  (4) 
Distant  flow  (l) 
Asymptotic  (5) 

The  results  show  that,  for  p  =  5  the  three  methods  give  almost  identi 
cal  answers.  For  p  <  5  the  asymptotic  solution  continues  to 


represent  the  exact  solution  well  for  large  even  for  p  as  small 
as  1.111.  The  distant  flow  approximation  (just  the  first  term  of  the 
asymptotic  series)  is  far  less  satisfactory  in  this  region  and  for 
p  =  1.111,  £  =  10.0  is  m  error  by  a  factor  of  almost  2.  This  region 
of  large  5  however  represents  the  part  of  the  cylinder  well  away  from 
tne  point  source,  where  the  flow  is  directed  principally  along  the 
axis  rather  than  across  it,  and  the  transverse  foree  is  very  small. 
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For  the  smaller  values  of  g,  closer  to  the  point  source,  the 
asymptotic  series  is  much  less  satisfactory  and,  surprisingly,  its 
first  tens  alone  (the  distant  flow  approximation) rep  resents  the 
integral  far  more  accurately.  This  is  nost  clearly  seen  in  the 
Figure  M3,  Although  only  expected  to  apply  for  (3  >  5»  even  for  £ 
as  small  as  1  .4  the  distant  flow  approximation  is  in  error  by  only 
12^  for  small  g.  Since  the  small  values  of  g  correspond  to  points 
near  the  point  source,  where  the  transverse  forces  are  largest,  the 
distant  flow  approximation  is  a  far  better  approximation  for 
practical  use  than  the  more  extended  asymptotic  result  [l5]f  and  can 
be  applied  down  to  very  short  ranges  {{3  «  1.4). 
b.  The  Force  Distribution  Near  the  Bow  and  Stern 

While  the  strip  theory  force  distribution  can  reasonably  be 
expected  to  apply  over  most  of  the  length  of  the  ship,  near  the  bow 
and  stern  it  might  be  expected  that  flow  around  the  ends  could 
relieve  the  fluid  pressures  and  so  alter  the  force  distributions  in 
these  areas.  For  the  overall  still  water  vibration  problem,  the 
inaccuracies  found  earlier  (Section  L)  for  the  strip  theory  forces 
did  not  make  much  difference  to  the  frequencies.  However,  for  a 
charge  close  to  either  the  bow  or  stem,  most  of  the  force  acting  on 
the  target  will  be  applied  in  this  error  prone  region  and  it  seems 
possible  that  small  errors  could  lead  to  large  differences  in  the 
forcing  function.  Since  it  can  allow  for  such  flow  the  3-D 
approximation  should  not  be  so  likely  to  cause  error.  The  accuracy 
of  both  the  strip  theory  and  the  3-D  approximation  are  investigated 
here. 

The  simplest  idealisation  to  use  in  this  case  is  an  ellipsoidal 
target,  with  a  point  source  to  represent  the  explosion  bubble.  As 
in  Section  L  the  ellipsoidal  co-ordinates  (£»  p,  w)  are  used  (they 


Ji  a r. 


—  i  ta  MiminIrTm 
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are  shown  in  Figire  R4).  They  are  related  to  cylindrical  co-ordi¬ 
nates  (r,  d,  x)  by 

r  =  k(£2  -  iV*  (1  -  ii2)* 
x  =  k  £p 

d=(o  where  k  is  a  constant.  [M6] 

In  terns  of  the  ellipsoidal  co-ordinates  the  general  solution  of 
Laplaccfe  equation  for  the  velocity  potential  <f>,  is  of  the  fora 
co  n 

<f>  =  2  £  (A  cos  aco  +  B  sin  a»)(C  P*  (£)  + 

n  n  n  n 

n=o  a=o 

“n  Qn  «»  *!  W  t«7l 

where  T®  (x),  P*  (x),  Q*  (x)  are  associated  Legendre  functions  of  the 
n  n  n 

first  and  second  kinds  and  are  defined  here  as 

T*  (x)  =  (-1)*  (i  -  x2)^*  ~  P (x)  -  1  <  x  <  l 

n  dx*  n 

PB  (x)  =  (x2  -  1  y*  £-  pn  (x)  x2  >  1  [M8] 

n  dx 

Q*  (x)  =  <-l)"(x2  -  1)^  Qn  (x)  x2  >  i 

n  dx*  n 

For  the  problem  of  a  point  source  of  strength  e(t)  at  the 
position  (£  =  £ ^ ,  p  =  fij »  »  =  o)  near  the  ellipsoid  £  =  CQ»  the 

potential  can  be  written  as 

=  $1  +  <f> 

where  ^  is  the  potential  of  the  point  source  alone  and  represents 

the  disturbance  potential  due  to  the  presence  of  the  rigid  ellipsoid. 
* 

4>  is  the  potential  needed  to  reduce  the  velocity  of  the  fluid  at  the 


ellipsoid  io  zero  in  the  direction  of  normals  to  the  ellipsoid. 
According  to  [H9]  the  potential  ,  in  the  form  [M7],  it  given  fay 


♦»  =  *  (2n 


*  »  [«»«>  P„  <V  Pn(p>  P„  <•*,>  * 


2  .E  f-£tsH|2  <  «>  P"  «,>  ^  <>*>  ^  <•*,>  «•  ~J 


for  values  of  C  >  and 


*rl  *  l2"*»  [«„  «,)  p.  <»  p„  <»>  pn  <•*,>  ♦ 
n=o 

2  j  <-'>*  {ferSHf  <  ‘V  pI «)  ^  w  i"  (m,) 

Ms  I 

cos  m»j  ;  C  <  {, 

# 

Since  <#>  *ust  be  continuous  and  finite  everywhere  outside  the  ellipse 
C  =  C0  its  expansion  aust  be  of  the  fora 


*  *  I  2  •»  [«.  K) 

n=o  ** 

cos  m»J  ;  £  >  £ 


(C)  Pn  <n)  ♦  ^  <£  (£)  ij  (H) 


where  the  coefficients  a  .  b _  are  to  be  chosen  to  aake 

n  lui 


“n  ■  -  to  *  0  on  5  *  «o- 


Since  on  the  ellipsoid  4-  «  4?«  this  condition  becomes 

on  0£ 
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so  that 


i,  b  Q*  '(C)  =  -  (2n  »  1).  2(-l)*  |f"  -  “(If 

n  na  n  o  in  +  aj. 


«S  «.>  C«o>  S  <“.> 


(from  the  orthogonality  of  the  functions  cos  aa  and  T*  (u)  over  the 

n 


ranges  (0,  2*),  (-1,  +l)  respectively), 
whence 

•  I 

W  p„  <•*,)  p„  w 


*  =  S  (2n  +  1 ) 

n=o 
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P JO-QJO  n  ° 
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n  <<«»>/ 


2 1  «!! «,)  ^  (m,)  %  w 


2 

v*  t (a  ■) l > 


a=1 


?;  to  -  q*  (o  x  «•  - 

n  oa 
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With  the  above  definitions  of  P  and  Q  the  Wronskian  relations  are 

n  n 


«> «;  *  «>  -  c  «> «>  -  {r^iH  p1:; 


so  that  on  the  ellipsoid  £  =  £  the  potential  reduces  to 


A  computer  routine  was  written  to  evaluate  this  series  at  a 
specified  set  of  values  of  x  along  the  ellipsoid  for  given  values  of 
a,  b,  c  and  x  .  The  Legendre  functions  are  readily  generated  fro* 
the  usual  recurrence  relation  except  for  the  Q's,  which  have  arguments 
greater  than  unity.  To  avoid  the  divergent  solution  in  this  case  it 
is  necessary  to  use  the  common  technique  of  starting  with  Q^  =  0, 

j  =  1  and  working  back.  This  gives  a  value  Q*  for  Qq  and 


multiplying  all  the  generated  Q's  by  -y -  gives  the  corrected 

Qo 


values.  Many  terms  in  the  series  are  necessary  as  the  charge 


approaches  the  ellipsoid  and  for  a  =  10,  b  =  1,  c  =  1.5,150  terms 
are  necessary  to  achieve  reasonable  accuracy  (O.ljQ,  For  a  =  10, 
b  =  1 ,  c  =  2.5,  50  terms  are  sufficient. 

The  results  given  fay  the  program  are  compared  in  Figure  15  with 
the  results  produced  by  both  strip  theory  and  fay  the  3-0  flow  tech¬ 
nique.  The  strip  theory  result  corresponding  to  [Ml 2]  is 


f  ,  •  (x) 

strip 


b2e  0  -  ■JJ2) 

2  -  -x>i3/2 


[M13] 


2xpe  Jc  +  (Xj  -  x)  y 

The  3-0  flow  result,  in  the  notation  of  Section  L  (equation  [L17])  is 


=  [b*  +  BA  u 

2xpe  L  -*  ^ 


[M14J 


^  22  • 

where  B  is  the  diagonal  matrix  (b^/l  )  and  uQ  is  given  by 


. 

u  . 
01 


Cc 
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Vac  result?  show  chat  in  fact  there  is  little  difference  between 
the  exact  result,  strip  theory,  and  the  3-D  approximation  when  the 
charge  is  at  one  end  of  the  ellipsoid,  even  when  it  is  very  close 
(c/b  =  1.5)*  Strip  theory,  although  slightly  less  accurate  than  the 
3-D  flow  approximation  near  the  peak  force  region  is  nevertheless 
quite  accurate.  The  table  shows  detailed  values  for  the  charge 
directly  opposite  the  end  of  the  ellipsoid.  The  two  3-D  flow  sets 
of  values  refer  to  the  results  given  by  equation  [ill  4]  above  firstly 
with  A  in  place  of  A  and  secondly  direct.  Although  significant  in 
the  vibration  case,  the  correction  is  everywhere  very  small  for  the 
explosion  induced  force  distribution. 

The  table  also  includes  values  for  the  case  of  the  charge 
opposite  the  centre  of  the  ellipsoid  (illustrated  in  Figure  M6). 

Again  strip  theory  is  a  very  good  approximation,  even  for  the  very 
close  geometry  but,  surprisingly,  in  this  case  the  3-D  approximation 
is  rather  less  satisfactory,  particularly  opposite  the  charge.  The 
difference  cannot  be  attributed  to  the  coarseness  of  the  dipole 
distribution  representation  since  doubling  the  number  of  constant 
dipole  segments  (giving  forty  segments  over  the  ellipsoid  length) 
confirmed  the  values  shown.  Considering  the  closeness  of  the  charge 
however,  even  the  3-D  flow  results  are  quite  reasonable. 

The  analysis  shows  that  as  regards  the  forcing  function  the 
strip  theory  and  the  3-D  flow  are  roughly  comparable,  strip  theory 
giving  better  results  near  the  target  centre  airi  3-D  flow  being 
slightly  better  near  the  target  ends.  In  either  case  both  methods 
give  very  reasonable  results,  particularly  since  it  is  not  the 
actual  values  of  force  per  unit  length  which  are  important  but. 
integrals  of  this  function  along  the  target  length.  It  is  clear 
from  the  figures  that  the  inaccuracies  in  such  integrals  will  be 
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even  less  than  in  the  values  of  force  per  unit  length 


It  may  be  noticed  from  Figures  M5  and  Jl6,  or  from  Figure  M3, 
that  even  when  the  charge  is  very  close  the  force  distribution  on 
the  target  is  spread  over  a  considerable  length  of  the  target.  This 
means  that  the  overall  forces  will  only  be  marginally  affected  by 
local  distortion  of  the  target  and  also  that  even  when  the  charge  is 
very  close  the  loading  is  never  really  equivalent  to  a  point  load. 

In  terns  of  the  normal  modes  excited, this  wide  load  region  means  that 
the  higher  modes  will  be  much  less  important  than  for  a  point  load, 
c.  Bernoulli  Pressures 

The  full  equation  for  the  pressure  in  a  fluid  in  unsteady  motion 
is 


d<J>  j,  2 

p  =  p_  +  P  at  “  3pu 

where  p  is  the  ambient  pressure,  u  the  fluid  velocity  and  <f>  the 


velocity  potential.  Throughout  the  theory  so  far  it  has  been  assumed 


that  the  inertial  pressure  p  -|^  was  much  more  important  than  the 

•Bernoulli'  pressure  !jpu  .  It  has  long  been  accepted  that  this  is 
certainly  a  reasonable  assumption  so  far  as  the  still  water  vibration 
modes  are  concerned  but  it  has  yet  to  be  established  for  the  flow 
around  the  bubble.  As  described  earlier  the  undisturbed  flow  around 
the  bubble  can  be  represented  quite  well  by  the  source  potential 

♦  -  ± 
o 

The  inertial  pressure  is  therefore 

•  • 

•  V 

pinertial  “  Pm>  “  4*r 


and  the  Bernoulli  pressure  is 


pBern  =  ^  =  * 


.2 


= 


2_4 


4xr  J  32-7C  r' 


The  ratio  of  the  two  pressures  is  therefore 


pBern  _  V2 

* *****  —  • 

pinert  8xrV 

On  using  the  bubble  scaling  factors  L  and  T  defined  in  Section  H, 
equation  [H5],  and  the  corresponding  non-dimensional  variables  x  and 
t,  this  becomes 


n  *  *2 

pbcrn  Lr_  (xJ) 

pinert  6r^  ^3 

Clearly  there  are  times  during  the  bubble  motion  when  jc*  is  almost 
»• 

zero  and  y?  is  large  so  that  the  inertial  term  dominates.  Equally 

*3  *3 

there  are  instants  when  x  is  zero,  and  at  these  x  is  obviously  a 
maximum  so  that  the  stagnation  term  dominates.  What  is  important  is 
the  relative  magnitudes  of  the  maxima  of  x^  and  y?  and  their 
durations. 

From  the  bubble  equations  given  in  Section  H,  the  volume 
acceleration  will  have  a  short  duration  but  very  high  maximum  at  the 
bubble  minimum,  and  a  long  duration  much  lower  secondary  'maximum' 

(a  suction)  at  the  bubble  maximum.  The  equations  give 


•  • 


max 


_2 
"  2 


k-zJl 

2 

x 

m 


I 
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so  that,  at  the  bubble  minimum 


x3  ~  -2  ".jiY— T  1 )—  ~  200  at  moderate  depths 
*  k“(y  -  1) 

and  at  the  bubble  maximum, 


-4  at  moderate  depths. 


For  must  be  zero  and,  from  equation  [Ht2]  this  will  occur 

max 


when  x3  »  <4.  Then,  from  equation  [H6], 

(x'*)  —  1 .75  again  at  moderate  depths, 

max 

On  using  these  values  it  is  clear  that 


pBern 

pinert 


,  a3 
l  _m 

7  r3 


where  a  is  the  maximum  bubble  radius 
m 


(normally  between  0.9L  and  0.94L) .  Hence,  the  maximum  Bernoulli 
stagnation  pressures  are  much  smellier  than  the  bubble  underpressures 
even  at  very  short  range  and  decrease  very  rapidly  compared  to  the 
inertial  pressures  at  all  greater  ranges.  In  addition,  the  duration 
of  the  Bernoulli  pressures,  being  related  to  the  period  when  the 
bubble  volume  expansion  rate  is  high,  will  be  much  less  than  that  of 
underpressure  phase  and  so  will  be  of  much  lesser  importance. 
Compared  to  the  peak  bubble  overpressures  the  stagnation  pressures 
are  completely  negligible.  The  stagnation  pressures  can 
consequently  be  ignored  at  all  ranges  down  at  least  to  the  range  at 
which  the  bubble  touches  the  target  at  its  maximum. 
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N  Modification  of  the  Forcing  Function  Due  to  Local  Damage 

The  discussion  in  Section  M  showed  that  the  forcing  (Unctions 
described  in  Section  H(  based  on  the  'distant  flow1  approximation, 
are  valid  down  to  very  short  ranges.  The  discussions  however  assumed 
throughout  that  the  target  remained  sensibly  rigid.  Since, 
particularly  at  short  ranges,  the  magnitudesof  possible  relative 
displacements  of  the  target  cross-sections  are  small  compared  to 
the  fluid  particle  displacements,  the  rigid  cross-section  assumption 
should  be  valid  unless  very  severe  damage  occurs.  The  type  of  severe 
damage  which  would  most  dramatically  change  the  fluid  flow  fields  is 
hull  rupture  since  the  presence  of  a  hole,  at  a  constant,  sensibly 
zero  pressure,  will  divert  part  of  the  flew  into  the  hole  instead  of 
out  around  the  target. 

An  approximation  to  the  local  flow  near  such  a  hole  is  obtained 
by  considering  the  flow  due  to  a  point  source  on  the  axis  of  a 
circular  hole  in  a  flat  plate.  This  problem  has  been  considered  by 
Haywood  [NlJ,  using  the  oblate  ellipsoidal  coordinates  (£,g,o>)  shown 
in  Figure  Nt •  These  coordinates  are  related  to  the  cylindrical 
coordinates  also  shown  by  the  equations 

2h  o  h 

r  =  a(l  +  C)  0  -  €  ) 
x  =  a£g 
6  -  « 

The  surfaces  £  =  const  form  a  family  of  confocal  oblate  ellipsoids 
and  the  surfaces  g  =  const  fona  the  orthogonal  family  of  confocal 
hyperboloids  of  one  sheet.  The  hole,  of  radius  a,  is  the  surface 
£  =  0.  The  plate  around  it  is  the  surface  g  =  0.  The  axis  through 
the  centre  of  the  hole  is  the  line  g  =  1  (a  degenerate  hyperboloid). 
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Haywood  found  that  the  velocity  potential  for  the  flow  produced 
by  a  source  of  strength  e  at  a  distance  xq  fro*  the  centre  of  the 

hole  (which  is  assuaed  to  be  a  surface  at  constant  zero  potential) 
can  be  written 


*  =  fj-  -  -  TZ  *F(C»  5)  +  FU’  (N,) 

where  R}  and  R2  are  the  distances  froa  the  source  and  its  iaage 
respectively.  The  function  F(£,  £)  is  given  by 


f«,  a  =  % 


rl  r  1 4 

J  (C  -  t)H  (Co  -  t?)2  +  (1  ♦  t2)(l  -  S2)  dt  (N2) 

i  *- 


the  source  being  at  the  point  (C0  =  *(/*•  5  =  0*  The  first  two 

terns  in  equation  (Nl)  give  the  usual  potential  for  a  source  near  a 
free  surface.  The  third  tera  represents  the  disturbance  in  the  flow 
due  to  presence  of  the  plate,  which  acts  as  a  baffle. 

The  integral  (N2)  can  be  written 


F(£.  5)  =  jsinh"1  (a,  -  ip,)  ♦  sinh"1  (-a,  -  ip,)] 
F<C,  O  =  -  2^  jsinh*1  (aj,  -  ip2)  +  sinh"1  (h^  -  ip2)J 


(C£0  ^  5)(C0  -  U) 

where  ai =  — r — 7* — pr  ; 
a.+no  ♦  c2)  (i  -  ? ) 

(c0  -  co2  ♦  (i  ♦  c2hi  -  s2) 

Pl  =  0  0  ^  ^ 
a  ♦  no  ♦  c2>  o  -  r) 
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(U0  -  S)U0  +&) 

2  2  ^  J2% 

0  ♦  no  +  o  0  -  c) 


<Cc  +  SC)2  +  (i  ♦  C2)(1  -  €2) 


0  +  nd  +  q  o  -  r) 


If  sinh”  (a  -  ip)  is  written  in  the  fora  y  +  i&*  then  y  and  5  are 
given  by 

2sinh2  y  =  |(i  -  a2  -  p2)  +  4a2J  -  (1  -  a2  -  p2) 

(K5) 

and  2  sin26  s  (1  +  a2  +  p2)2  -  J(l  +  a2  +  p2)  -  4p^j 

where  positive  values  are  to  be  taken  for  the  square  roots.  Clearly 
both  y  and  6  are  multi-valued;  both  being  indeterminate  in  sign  and 
6  also  being  indeterminate  to  a  multiple  of  x.  The  correct  values  to 
be  used  must  be  selected  to  ensure  that  <f>  is  continuous  within  the 
fluid.  From  (N5),  y  changes  sign  only  when  a  =  0  so  that  the  real 
parts  of  sinh  1  (a  -  ip)  and  sinh  1  (-a  -ip)  are  equal  and  opposite 
and  F(C,  £)  and  F(£,  -£)  are  both  purely  imaginary.  The  potential  <f> 
is  therefore  real  as  required.  Similarly,  6  can  only  change 
sign  when  p  =  0  and,  since  p  is  always  >  0,  cannot  change  sign. 

It  is  therefore  restricted  to  the  range  0  <  8  <  x.  This  leaves  open 
which  of  the  ranges  (0,  %/2 )  and  (x/2,  x)  8  occupies.  It  is  clear 
from  (N5)  that  8  can  change  from  one  of  these  ranges  to  the  other 
only  when  a  =  0  and  then  only  if  6  >1.  For  F(£,  g),  =0  when 


g  »  and  then  p}  = 


k2  - 

and  so  is  always  <  1 .  8,  will 

««  * 


therefore  be  confined  to  the  range  (0,  r/2) .  For  F(C»  -£)»  a2  =  0 

2  ^ 

<«o  +  O 

when  S  =  CC  and  then  ■  ■  which  is  always  >  1 . 

°  “  o 

0  -  c> 

Consequently  &2  changes  from  one  of  the  ranges  (0,  x/2)  and  ( r/2 ,  it) 

to  the  other  when  a2  changes  sign.  Since  <f>  must  be  zero  for  £  =  0, 

the  range  (0,  %/2 )  corresponds  to  a  <  0. 

These  results  may  be  summarised  by 


where  6^  and  &2  are  the  positive  angles  given  by  equations  fN4)  and 
(N5),  with  5^  in  the  range  (0,  t/2 )  and  &2  in  the  range  (0,  r/2)  or 
(*/2t  x)  according  as  is  less  than  or  greater  than  zero.  A 


simple  verification  is  provided  by  two  special  cases: - 

(i)  a  -4  oo  (  so  that  there  is  no  plate,  only  a  source  near  a 
free  surface). 

Then,  for  moderate  values  of  xq,  x  and  r 


_  ,  ,  ->  (x°  -  *  *  -* 
1  r  J  ar 


and  so  6j  0.  Also, 


_/X  »  ,x  \  2  2 

v  _  sC°-+  ■*)  .  A  ^  Lft  JLMi  XJ £■  o 


2  r  ’  r2 

so  that  &2  tends  to  either  0  or*,  Since  <  0,  &2  0 


R1  R2 


and  <f> 


which  is  the  usual  potential  for  a  source  near  a  free 


surface. 

(ii)  a  0  (so  that  there  is  no  hole  in  the  plate). 


Then 


a» 


Pt  0 


*2 


x 

o  +  X 

r 


P2  -4  0 


and  again  6^  0.  This  time  since  a,,  >  0,  62  %  and 

the  potential  <f>  has  the  limiting  value 


which  again  is  the  correct  potential  for  a  source  near  a 
rigid  plate. 

Neglecting  the  Bernoulli  pressures  (in  keeping  with  the 
previous  approximations)  the  pressure  on  the  plate  is  given  by 


P 


(N7) 


since  >  0  for  g  =  0. 

If  there  were  no  hole  in  the  plate,  5^  would  be  zero  and  the 
pressure  on  the  plate  would  be  simply 


P  = 


-  2££ 


'1 


The  effect  of  the  hole  is  consequently  to  reduce  the  pressures  by 
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the  factor  (t  -  — — )  from  the  value  for  the  plate  without  a  hole. 


Figure  N2  shows  the  values  of  p(r)/2pe  for  several  values  of  xy'a. 
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For  large  values  of  r,  the  factor  (t  -  — — rapidly  approaches  its 

A 

asymptotic  value 

R  =  1  -  7  tan”*  (^-)  (N8) 

•o  7v  A, 

o 

so  that,  except  very  close  to  the  hole,  the  pressures  everywhere  are 
given  by 

P  =  Po  R» 

where  pQ  is  the  pressure  ignoring  the  hole.  This  approximation 

provides  a  very  reasonable  method  for  allowing  for  the  effect  of  the 
hole.  All  forces  on  the  target  are  calculated  as  though  the  forcing 

function  were  reduced  in  strength  by  the  factor  R  . 

00 

The  factor  R  has  a  simple  physical  interpretation.  The  amount 

OO 

of  fluid  flowing  through  the  hole  per  unit  time  may  be  shown  by 

integration  across  the  hole  to  be  given  by  8  eptan”*  (a/xQ).  Since 

the  total  mass  of  fluid  flowing  from  the  source  per  unit  time  is 
43tpe  the  fraction  of  this  which  enters  the  hole  is 


and  the  remaining  fraction  of  fluid,  which  avoids  the  hole,  is 
simply  R  ,  The  source  and  hole  together  consequently  appear,  at  a 

CO 


distance,  as  a  source  of  reduced  strength  eR  .  Figure  N3  shows  the 

00 


variation  of  R  with  the  ratio  x /a 

00  Q 


This  approximation  should  be  quite  reasonable  when  the  charge  is 


close  to  the  target  but  as  the  charge  standoff  is  increased  the 
curvature  of  the  cylinder  will  begin  to  affect  the  quantity  of  water 
which  is  deflected  through  the  hole.  It  will  become  easier  for  the 
fluid  to  flow  around  the  cylinder.  In  consequence  less  water  will  be 
lost  through  the  hole  and  the  factor  R  will  become  slightly  larger 
than  for  the  corresponding  standoff  from  a  flat  plate.  The  magnitude 
of  this  effect  is  readily  estimated.  The  new  value  of  R  must 
obviously  be  less  than  the  value  which  would  pertain  if  there  were  no 
plate  at  all  to  deflect  flow  into  the  hole,  just  a  hole  in  the  fluid. 
This  extremely  artificial  case  can  be  solved  fairly  easily. 

Bryant  [N2]  determined  the  potential  due  to  a  point  source  on  the 
axis  of  a  rigid  circular  disc  to  be 


on  the  surface  of  the  disc  facing  the  source.  In  the  problem  to 
be  solved,  the  same  boundary  condition  over  the  rear  of  the  disc 

applies,  ie  ^  =  Otbut  over  the  front  face  the  condition  should  now 

be  <t>  =  0.  The  potential  for  this  problem  may  be  found  from  Bryant's 
by  adding  a  'disturbance'  potential  <f>2  given  by 


<t>2  =  e  S  an  Qn  (i{)  Pn  U)  (N10) 

n=o 

where  the  coefficients  aR  are  chosen  to  satisfy 
<f>2  =  “  ^or  front  face  of  the  disc 
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and 


(Nil) 


&<t>2 

■5 —  =0  on  the  rear  face, 

on 

This  nixed  type  of  boundary  condition  is  difficult  to  satisfy 
exactly  but  if  the  series  for  <f>2  is  assumed  to  be  convergent  then  the 

first  few  coefficients  can  be  found  approximately  by  simply 
satisfying  the  appropriate  boundary  conditions  at  an  arbitrary  set 
of  collocation  points.  From  the  nature  of  the  problem  it  is  likely 
that  these  collocation  points  would  be  best  chosen  as  the  zero's  of 
an  appropriate  Tschebyscheff  polynomial.  However,  for  simplicity  a 
set  of  equally  spaced  points  was  in  fact  used.  The  condition  on  the 
front  face  was  satisfied  at  n  internal  points  and  also  at  the  rim  of 
the  hole.  The  condition  on  the  rear  face  was  satisfied  at  n  internal 
points.  A  short  computer  program  was  used  to  solve  for  the  first 
(2n  +1)  coefficients  in  the  series.  The  coefficients  were  determined 
for  several  values  of  n  and  their  values  were  found  to  be  rather 
slowly  convergent  but  the  value  of  the  first  coefficient,  aQ  had 

settled  to  under  0.25/J  variation  for  n  =  10,  ie  for  a  21  term 
approximation  to  the  solution. 

With  the  coefficients  in  (N10)  known  the  quantity  of  fluid 
flowing  into  the  hole  is  most  easily  found  by  considering  the  flow 
at  very  large  distances  from  the  origin.  Then  the  ellipsoids 
C  =  const  are  very  nearly  spheres  of  radius  a£  so  that  the  outward 
flow  across  such  a  surface  per  unit  time,  due  to  the  disturbance 
potential  is  approximately 

2  2  1  ^$2  m  2  / 

4*  (a£)  =  -  4*  (aO  .  -  =  -  4*aei  Z  anC  Q„( iOP„( 5) . 

n=o 
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1 .3.3 


Nov*,  for  large  £,  Q  (i£) 


J'i  ,-(n-l)  „  ,,-(n+3h 

—  (T„  rn f  * 0  (?  > 


(referenceN3,  art  106).  so  that  in  the  limit  as  £  — >  »,  C  Qn  "■>  0 

for  all  n  *  1  and  the  outflow  at  »  is  siaply 
4*aei  a 


The  outflow  due  to  the  original  potential  ^  is  the  source  outflow 


4«e  so  that  -4r.aeiaQ  represents  the  decrease  in  outflow  due  to 

the  hole.  The  reduction  coefficient  in  this  case  is  then  simply 
R  =  l  ♦  iaa 

o 

This  function  is  compared  with  R  in  Figure  N3.  For  very  small 

00 

values  of  x^a  the  source  sees  only  the  hole,  not  the  plate  or 

fluid  surrounding  it  and  so  R  and  R  are  identical.  As  x 

«  o 

increases  the  effect  of  the  difference  becomes  apparent  but  is 

always  relatively  small.  Since  R  is  an  upper  limit  for  the  true 

reduction  coefficient,  which  would  be  expected  to  be  much  closer 

to  R  than  to  R,  use  of  R  for  cylinders  seems  completely 
00  00 


justified. 
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GEOMETRY  FOR  A  SOURCE  NEAR  A  HOLE  (N  A  PLATE 


ASYMPTOTIC  REDUCTION  FACTOR  AS  A  FUNCTION 
OF  STANDOFF  /  HOLE  RADIUS  RATIO 

FIGURE  N  3 
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0  The  Effect  of  Bubble  Migration 

Throughout  the  entire  discussion  so  far  on  the  limitations  of 
the  theory  it  has  been  assumed  that  the  bubble  remained  fixed  in 
space  throughout  its  pulsations.  For  small  charges  of  the  type  used 
in  Chertock's  experiments  this  is  a  reasonable  assumption  since, 
except  in  special  circumstances  (eg  an  underpressure  tank),  the 
pressure  difference  across  the  bubble  is  small  compared  to  the 
absolute  hydrostatic  pressure.  For  larger  charges,  of  the  sizes 
likely  to  affect  ships,  the  pressure  differences  will  be  very  much 
larger  but  the  hydrostatic  pressure  will  not  change  greatly.  For 
example,  consider  the  model  in  Chertock's  experiments.  It  was 
10  ft  long  and  the  maximum  bubble  radius  was  about  6  ins  when  the 
charge  was  suspended  l6  ins  beneath  the  water  surface.  On  full 
scale  the  target  might  be  40  times  as  big  so  that  a  scaled  bubble 
would  have  a  20  ft  radius  and  its  centre  would  be  50  ft  beneath  the 
surface.  In  the  model  case  the  ratio  of  (pressure  difference 
between  the  bubble  top  and  bottom)  over  (pressure  at  bubble  centre) 
is 

*^)  .  ,  =  vj  ,1'*^  "f*T  vr  =  0,031  where  atmospheric  pressure 

p  moM  33  rt*  1.3  ft  is  equivalent  to  .  33  ft 

head  of  water 

On  full  scale  the  comparable  ratio  is 
■p^full  scale  a  33*+  50  =  0,48 

Clearly  in  this  case  the  relative  buoyancy  forces  acting  on  the 
bubble  will  be  sixteen  times  greater  than  on  the  small  scale  and 
migration  will  be  more  important. 
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The  mechanism  of  bubble  migration  is  ver;-  well  understood 
(reference  01 ).  During  the  early  expansion  very  little  migration 
occurs  since  the  bubble  is  snail  and  the  hydrostatic  force  on  it 
(its  displacement)  is  similarly  small.  During  its  underpressure 
phase,  when  it  is  large,  the  buoyancy  force  is  very  large  but  so  is 
the  effective  mass  of  the  bubble;  although  the  gas  within  the  bubble 
has  a  negligible  mass  compared  with  the  mass  of  water  displaced  by 
it,  the  added  mass  associated  with  the  flow  of  water  around  the 
bubble  is  not  negligible  -  for  a  spherical  bubble  it  is  half  the 
mass  of  the  displaced  water.  In  oonsecpence  the  bubble  accelerates 
upwards  but  does  not  have  time  to  achieve  a  significant  displace¬ 
ment  or  velocity  before  the  bubble  contraction  occurs.  It  does 
however  have  time  to  acquire  a  considerable  momentum  (actually 
contained  in  the  water  around  it).  During  the  contraction  phase  the 
mass  of  water  displaced  is  reduced  to  a  very  small  fraction  of 
maximum  displacement  and  the  added  mass  is  naturally  similarly 
reduced.  Since  the  momentum  acquired  during  the  underpressure  phase 
is  still  associated  with  the  bubble  the  reduction  in  added  mass 
implies  an  equivalent  increase  in  vertical  velocity  and  the  bubble 
shoots  upwards  when  near  its  minimum.  During  this  phase  it  distorts 
greatly,  the  bottom  surface  iiwerts  and  forms  a  jet  which  penetrates 
through  the  bubble,  producing  a  toroidal  shape.  During  re-expansion 
the  bubble  regains  an  approximately  spherical  shape  and  its  velocity 
falls  as  its  added  mass  rises.  Thus  the  bubble  is  almost  stationary 
during  the  phase  when  the  vertical  forces  on  it  are  large,  but 
migrates  rapidly  later  when  the  vertical  forces  are  least.  During 
the  rapid  migration  phase,  since  the  vertical  force  is  small  the 
momentum  is  almost  constant  so  that  the  product  of  bubble  volume  and 
vertical  velocity  is  almost  constant. 


Trie  wain  and  most  obvious  effect  of  migration  is  to  bring  the 
point  of  emission  of  the  first  bubble  pulse  much  nearer  to  the  surface 
of  the  water.  Since  this  pulse,  as  discussed  earlier,  can  be  almost 
twice  as  effective  as  the  initial  expansion,  this  can  considerably 
affect  the  motion  of  a  target  above  the  bubble.  There  are  however 
other  effects.  The  peak  pressure  during  the  bubble  pulse  is  well 
knewn  to  be  considerably  reduced  by  vigorous  migration  (reference  Ol). 
The  reduction  in  pressure  however  is  accompanied  by  an  increase  in 
duration  so  that  the  impulse  is  less  affected.  In  addition,  the 
motion  of  the  bubble  causes  a  dipole  pressure  field  around  it  in 
addition  to  the  usual  radial  field.  This  dipole  field  was  first 
investigated  by  Taylor  [02]  who  showed  how  very  intense  it  can  be, 
although  of  course  it  decays  much  more  rapidly  with  distance  from 
the  bubble  than  the  radial  field  and  is  seldom  observed  in  pressure 
gauge  measurements. 

Another  effect  due  to  gravity  is  distortion  of  the  bubble  from 
its  spherical  form.  Such  deviations  have  been  the  subject  of  a 
number  of  studies,  eg  [03]  and  [04],  but  in  the  main  these  have  been 
based  on  perturbation  techniques  and  have  only  applied  to  the  early 
bubble  deformations,  Quito  recently,  purely  numerical  techniques 
have  been  used  to  describe  the  flow  around  pulsating  bubbles, 
particularly  around  collapsing  cavitation  bubbles.  These  have  shown 
clearly  [05]  the  sort  of  distortions  suffered  by  a  bubble  collapsing 
in  an  external  pressure  gradient,  although  in  the  case  of  [053  the 
gradient  was  due  to  the  presence  of  a  rigid  wall  rather  than  a 
gravity  field.  These  distortions  are  related  to  multipole  terms  in 
the  flow  field  of  higher  order  than  dipole.  These  higher  order  terras 
decay  with  distance  even  more  rapidly  than  the  dipole  terms  and  90 
are  neglected  here.  The  distortions  lead  to  drastic  changes  in  the 
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motion  of  the  bubble  during  the  collapse  phase  when  near  the  minimum 
and  a  significant  redistribution  of  energy  occurs  at  this  time. 

Since  such  effects  cannot  at  present  be  treated  analytically,  they 
are  allowed  for  here  only  on  an  empirical  basis. 

With  the  geometry  of  Figure  01,  assuming  first  that  the  explosion 
is  deep  so  that  the  effect  of  the  water  surface  can  be  neglected,  the 
velocity  potential  for  the  migrating  bubble  will  be  approximately 


“2  C0S  ®1 
r1 


(01) 


While  the  bubble  is  spherical  this  will  represent  the  motion  exactly 
but  late  during  the  bubble  collapse  phase,  when  it  begins  to  distort, 
higher  order  terms  will  also  be  present.  Since  these  will  affect  the 
flow  only  very  close  to  the  bubble  surface  they  are  neglected  and 
equation  (01 )  is  assumed  to  apply  throughout  the  motion,  the  bubble 
being  assumed  to  remain  spherical  at  all  times.  The  strengths  e^  and 

e 2  are  then  readily  found  to  be  related  to  the  bubble  volume  V(t)  and 
velocity  v(t)  by 


t}  =  V/4 %  and  e2  =  (02) 

From  equation  (01 )  the  vertical  velocity  due  to  the  pulsating 
migrating  bubble  is 


y  *  Fry  '  r  3 

i  i 


r2 
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(03) 


Since  the  bubble  is  rising  with  velocity  v,  the  vertical  accelera¬ 


tion  Uy  is  given  by 


.  duv  \  «,  (y +  °)  \ 

Uy  =  '  v*sr =  n - 
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(04) 


The  nature  of  the  four  terms  is  easily  seen  by  substituting  the 

expressions  (02)  for  e1  and  e2*  The  first  term,  proportional  to 

•  •• 
is  therefore  proportional  to  the  bubble  volume  acceleration,  V, 

and  so  in  principle  is  the  same  as  for  the  non-migrating  bubble 

•  • 

(although  the  value  of  V  will  be  changed).  The  second  term  is 


proportional  to-^(Vv).  Now  the  mass  of  the  bubble  itself  is 

negligible  compared  to  the  added  mass  of  the  water  moving  around  it, 
which  for  a  spherical  bubble  is  just  one  half  the  displaced  water 
mass.  The  quantity  Vv  is  therefore  proportional  to  the  momentum 
associated  with  the  moving  bubble  and  the  term  e^  is  consequently 

proportional  to  the  rate  of  change  of  bubble  momentum,  and  so  also 
to  the  force  acting  on  the  bubble.  This  force  is  simply  the 
buoyancy  force  due  to  the  gravity  induced  pressure  gradient  in  the 
water,  ie  pVg.  The  second  term  is  therefore  proportional  to  the 
bubble  volume.  It  will  be  negligible  in  the  early  expansion  and 
again  later  near  the  bubble  minimum.  It  will  be  greatest  during 

the  long  expansion  phase.  The  third  term  is  proportional  to 

•  • 
e^  v  <*  Vv  and  will  therefore  be  greatest  when  V  and  v  are  both 
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large,  ie  late  in  the  collapse  phase  and  again  early  in  the 
subsequent  re-expansion.  It  will  be  very  small  during  the  under- 
pressure  phase  since  then  both  V  and  v  are  very  saall,  and  will 
also  be  small  at  the  bubble  minimum  since  V  is  zero  then.  The 

fourth  term  is  proportional  to  e0v  <*  Vv  <x  a  /V  where  a  is  the 

vertical  aoaentum  of  the  bubble.  It  will  consequently  be  largest 
at  the  bubble  ainiaua. 

On  assuming  that  the  bubble  remains  spherical,  the  kinetic 
energy  equation  given  in  section  H,  when  applied  to  the  moving 
bubble,  becomes 


E  =  2xpa3a2  +  4xpa3gZ  +  y  . -  * -  + 

a*-’) 

~  xpa3v2  (05a) 

Equating  the  buoyancy  force  to  the  rate  of  change  of  aoaentum 
gives  the  migration  equation 


|  *pa3vj  =  j  *pa3g  (05b) 

In  terms  of  the  non-diaensional  variables  x  and  t  defined  in 
section  H,  and  with  Z/l  =  £,  these  equations  become 

i2  *  £2/6  ♦{/£„-  |i  -  k/x3(Y"1  ’j/x3 


•Sj  (X3  C)  *  -  (06) 

and  the  four  terms  in  the  expression  (04)  for  the  fluid  accelera¬ 
tion  may  be  re-written  (for  the  surface  y  =  0) 
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In  each  of  these  terns  the  acceleration  is  aeasured  in  units  of  g, 
the  acceleration  of  gravity.  The  expressions  (07)  have  all  been 
multiplied  by  the  factor  2  to  allow  for  usual  doubling  effect  of 
the  free  surface.  Since,  for  rapid  notions  where  surface  waves  can 
be  neglected,  the  boundary  condition  to  be  satisfied  is  <f>  s  0,  the 
condition  is  satisfied  by  a  negative  inage  above  the  surface  which 
produces  the  assumed  enhanced  acceleration. 

The  Free  Surface  Effect 

The  free  surface  is  well  known  to  affect  not  only  the 
accelerations  near  it  but  also  the  notion  of  the  bubble  itself  if 
it  is  within  a  few  bubble  naxiaua  radii.  The  velocity  potential  in 
such  cases  is  given  approxinately  by 


cos 


(OS) 


for  the  geometry  of  figure  01 .  In  this  case,  as  Taylor  showed,  the 
source  and  dipole  strengths  are,  to  order  (aVd^), 
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The  kinetic  energy  of  the  flow,  to  the  same  order,  is 
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and  the  energy  equation  becomes 

E„  =  2*paV  (l  -  fj)  *  j  f»3Z2  ♦  f  p«3ii  4  *  f'P"3*1 


«a^p 
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where  p  is  the  pressure  within  the  bubble. 

The  migration  equation  may  in  this  case  be  derived  from  the 
Lagrangian  form  of  the  equations  of  motion,  ie 


(T  -  V)] 
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In  terms  of  the  previous  non-dimensional  variables  the  resulting 
equations  are 
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Since  Z  =  d  +  h,  h  being  the  static  head  of  water  equivalent  to 
atmospheric  pressure,  6  =  £  -  where  =  h/L.  For  computational 

purposes  these  equations  are  best  reduced  to  the  equivalent  set  of 
first  order  equations 
x  =  a 


X  =  -3apr-+  “  +  (3cr2  +  Xff)j  (01 1 ) 

where  two  parameters  a  and  (3  have  been  introduced  to  enable  the 
same  equations  to  cover  three  distinct  cases  of  idealised  bubble 
motion*  The  appropriate  values  for  a  and  p  are 


The  appropriate  initial  conditions  are  those  used  earlier  in  the 
non-migrating  case,  ie 


<7=0 

X  *  0 


close  to  unity  in  value,  is  vital  to  the  successful  integration  of 
the  equations.  A  computer  program  to  integrate  these  equations 
has  been  written,  using  a  standard  Runge-Kutta  integration  routine 
to  integrate  the  first  order  differential  equations  (Oil).  Because 
the  variables  change  very  rapidly  near  the  bubble  minima,  but  are 
relatively  static  during  the  underpressure  phase,  a  constant  time 
increment  in  the  integration  routine  is  very  inefficient  and  the 
program  has  been  written  to  select  its  own  step  length. 

Corresponding  to  the  use  of  equations  (Oil)  for  the  bubble 
behaviour,  the  four  component  accelerations  at  the  free  surface 
become 


(012) 


Figure  02  shows  the  bubble  motion  and  the  surface  accelerations 
above  the  bubble  for  a  typical  large  explosion  bubble.  The  figure 
compares  the  non-migrating  results  with  those  given  by  equat ions  (06) 
and  (07)  for  the  migrating  case  and  v/ith  the  results  from  (Oil)  and 
(01 2)  for  the  same  case  with  surface  effects  included.  In  the 
figure  only  the  results  prior  to  the  bubble  minimum  are  given  for 


the  migrating  cases  since  the  bubble  equations  break  down  at  the 
minimum  due  to  proximity  to  the  free  surface.  The  necessity  of 
allowing  both  for  migration  and  for  the  free  surface  is  clear. 
Inclusion  of  a  Drag  Term 

In  all  the  above  equations  the  accelerations  depend  directly 
on  the  migration  velocities  of  the  bubble,  as  deduced  from  ideal 
flow  around  a  spherical  bubble.  It  is  well  known  however,  eg 
[o6],  that  such  migrations  are  considerably  larger  than  the 
measured  values  and  the  surface  acceleration  deduced  from  them  will 
be  correspondingly  in  error.  Since  it  is  impractical  to  attempt 
to  describe  the  full,  distorting  flow  around  the  bubble  except  by 
very  extensive  numerical  methods  it  is  necessary  to  introduce  an 
empirical  correction  to  the  migration  equations. 

The  most  obvious  such  correction  would  be  to  give  the 
parameter  a  introduced  into  equation  (Oil)  a  suitable  value  less 
than  unity.  Unfortunately  this  simple  technique,  through  equations 
(012)  would  affect  the  f£  component  of  the  acceleration.  This 
component  is  a  maximum  when  the  bubble  itself  is  a  maximum  and  is 
scarcely  moving,  ie  when  the  ideal  flow  conditions  are  most  nearly 
satisfied  and  it  would  be  expected  that  the  component  f2  should  be 

most  accurate.  Also,  it  would  affect  small  weakly  migrating 
bubbles  relatively  as  much  as  large  strongly  migrating  ones  since 
a  reduction  in  the  value  of  a  is  equivalent  to  a  reduction  in  the 
•  buoyancy  force  acting  on  the  bubble.  What  is  required  is  an  effect 
which  is  principally  ajjparent  when  the  bubble  is  migrating  most 
rapidly.  This  suggests  a  turbulent  wake  type  of  drag  force  of  the 
form 

F  =  CD  (^  pv2  A) 


loO 


where  A  is  the  bubble  cross-sectional  area  transverse  to  the  raoKon 


and  Cp  is  a  drag  coefficient.  If  the  bubble  were  .  solid  sphere 

such  a  drag  term  would  certainly  exist  as  a  typical  bubble  from  a 
large  explosion  could  easily  be  5  ft  in  radius  when  near  its 
minimum  and  moving  at  up  to  200  ft/sec.  The  Reynolds  number  for 

g 

such  a  motion  is  of  the  order  of  10  »  Since  the  actual  bubble  is 
a  surface  of  constant  pressure,  the  mechanism  of  wa^e  formation  is 
not  so  obvious  as  for  a  solid  body  but  some  such  force  seems  likely. 
A  correct  physical  picture  is  not  however  required  of  an  eopirical 
correction.  It  is  sufficient  that  the  assumed  force  should  lead  to 
realistic  results. 

Assuming  the  existence  of  a  force  of  the  above  type,  the  motion 
of  the  bubble  will  be  expending  energy  at  the  rate 

dE  -  X  .  2  3  x  „  2  *3 

'  dt  =  Fv  =  2  PCD  "  ^  “  -2  P°D*  2 

Similarly  momentum  is  being  lost  to  the  wake  at  the  rate  F 
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Although  the  Lagrangian  derivation  of  the  bubble  motion  equations 
is  no  longer  available,  the  simpler  type  of  ^derivation*  due  to 
Taylor  leads  to  the  equations 
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These  equations  are  equivalent  to  the  first  order  equations 
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The  corresponding  component  accelerations  at  the  free  surface 
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The  computer  program  for  the  bubble  motion,  which  was  described 
earlier,  in  fact  integrates  equations  (013)  (all  the  previous  cases 
can  be  obtained  from  (01 3)  using  suitable  values  for  a,  {3  and  C^). 

A  series  of  calculations  were  carried  out  for  a  500  lb  charge  of  TNT 

at  a  depth  of  150  ft  using  a  series  of  values  for  C^.  The  results 

are  summarised  in  the  Table. 

Migration  Variables  for  Different  Drag  Coefficients 
for  a  500  lb  Charge  at  150  Ft 


Drag  Coefficient 

CD 

Migration  to 
First  Pulse 
(ft) 

P1 

Edissipated 

E 

0 

0.00 

16.9 

0.47 

0.00 

0.50 

15.0 

0.37 

o.»9 

1  .00 

15.3 

0.40 

0.22 

1.50 

12.3 

0.46 

0.19 

2.00 

11.3 

0.52 

0.15 

2.25 

10.9 

0.56 

0.12 

Experimental  Result 

11 .0 

0.63 

0.15 

In  the  table  p^  is  any  pressure  which  would  be  measured  at  large 

standoff  and  p^  ^  is  the  pressure  which  would  be  measured  at  the 

same  point  if  the  bubble  did  not  migrate.  In  the  near  field  the 
pressures  will  be  affected  by  the  dipole  field  but  at  large 
distances  only  the  radial  field  will  be  felt  and  the  ratio 


'"ToT  tends  quite  rapidly  to  a  limiting  value.  The  experimental 
P1 

value  for  p^  is  subject  to  considerable  scatter.  The  ’experimental’ 

value  for  p  ^  is  estimated  from  the  empirical  formula  for  bubble 
pulse  pressures.  The  table  shows  that  for  this  particular  charge 


1b3 


size  and  depth  the  value  =  2.25  makes  both  the  migration  and  the 

P, 

ratio  ■— approximately  correct. 

P1 

The  large  value  of  Cp  precludes  its  direct  physical  identifi¬ 
cation  with  a  normal  wake-type  drag  coefficient  since  its  value 
would  normally  be  fairly  close  to  unity.  However,  for  the  strongly 
migrating  bubbles  for  which  the  value  was  derived  the  bubble  would 
have  flattened  considerably  near  the  minimum  (as  is  observed  on 
small  scale)  so  that  the  area  to  which  is  related  would  actually 

have  been  rather  larger  than  assumed.  Had  a  nu.re  realistic  value 
for  the  area  been  used  a  lower  more  realistic  value  of  would 

have  resulted.  However,  a  complete  physical  picture  is  not  required 
here,  just  a  simple  mechanism  for  producing  the  correct  bubble 
migrations. 

One  of  the  equations  from  which  the  bubble  motion  equations 
were  derived  was  the  energy  balance  equation.  The  addition  of  the 
drag  terra  implies  that  the  bubble  loses  energy  continually  through¬ 
out  the  motion.  Most  of  this  energy  loss  naturally  occurs  at  the 
bubble  minimum  when  the  velocity  is  very  high.  It  is  well  known 
that  the  period  of  a  pulsating  bubble  is  quite  a  good  measure  of 
the  bubble  energy  and  this  fact  has  long  been  used  to  measure  the 
bubble  energy.  By  comparing  the  periods  of  the  first  and  second 
bubble  pulsations  estimates  have  been  made  of  the  total  energy  loss 
which  occurs  at  the  bubble  minimum  (reference  [06]).  The  energy 
still  retained  by  the  bubble  at  its  minimum  was  lculated  for  the 
coses  given  in  the  table  and  the  energy  dissipated  up  to  this  time, 
as  a  fraction  of  the  original  energy,  has  been  added  to  the  table. 
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The  experimental  value  for  the  energy  los-  for  this  charge  size  and 
depth  is  40%  of  the  original  bubble  energy,  of  which  1C#  is  estimated 
to  be  due  to  acoustic  radiation  (neglected  in  this  analysis).  -  If  it 
is  assumed  that  roughly  half  the  remaining  dissipation  occurs  before 
the  minimum,  this  gives  the  155»»  quoted  in  the  table,  for  the 
experimental  dissipation.  The  value  Cp  =  2.25  therefore  gives  about 

the  right  total  energy  dissipation.  It  is  of  interest  -hat  the  value 
Cp  =  1.00  corresponds  to  a  maximum  in  the  energy  dissipated. 

One  significant  point  regarding  the  energy  loss  is  that  although 
for  very  weakly  migrating  bubbles  an  energy  dissipation  mechanism  is 
quil i  clear,  the  same  is  not  true  of  strongly  migrating  bubbles.  For 
small  bubbles,  or  very  deep  ones,  which  do  not  migrate  significantly, 
the  bubble  surface  remains  fairly  spherical  until  late  in  the  collapse 
phase  when  the  bubble  pressure  rises  above  the  ambient  level  and 
starts  to  decelerate  the  bubble  surface.  The  acceleration  then  being 
directed  from  the  less  dense  f*uid  to  the  more  dense  one,  causes  the 
usual  Taylor  surface  instability  to  occur.  High  speed  films  show 
thc.t  the  instability  occurs  as  a  very  large  number  of  needle  like 
jets  being  projected  into  the  bubble.  Due  to  its  adiabatic  compres¬ 
sion  the  gas  will  be  very  hot  at  this  time  and  such  an  intense  spray 
of  cold  water  will  remove  much  energy.  For  strongly  migrating 
bubbles  this  f^rra  of  instability  does  not  seem  to  occur.  An  over¬ 
all  distortion,  first  into  an  oblate  spheroid,  occurs  instead  and 
later  the  bottom  inverts  and  penetrates  through  the  bubble  turning 
it  into  a  toroid.  Since  intimate  mixing  of  the  gas  and  water  no 
longer  occurs  an  alternative  mechanism  is  required  to  account  for 
the  observation  that  migrating  bubbles  dissipate  as  much  energy  as 
non-migrating  ones.  The  'wake*  theory  provides  such  a  dissipation 


mechanism.  The  energy  is  not  so  much  dissipated  as  transferred  to 
water  which  loses  its  association  with  the  bubble.  The  mechanism 
assumes  that  both  momentum  and  energy  are  transferred  to  water 
behind  the  bubble,  principally  when  the  bubble  is  small  and  its 
velocity  great.  As  the  bubble  begins  to  re-expand  and  slow  down, 
it  has  no  way  of  stopping  the  narrow  following  wake  which  could 
penetrate  it  and  form  a  vertical  water  jet.  The  large  plumes  which 
are  so  evident  on  the  surface  above  an  explosion  are  known  to 
originate  at  the  bubble  pulses  and  are  perhaps  due  to  such  jets. 

One  important  result  of  the  transfer  of  energy  and  momentum 
to  the  'wake'  of  the  bubble  is  that  the  bubble  motion  is  no  longer 
symmetrical  -»bout  the  minimum.  Figure  03  shows  the  displacement 
time  history  calculated  for  the  bubble  centre  for  the  500  lb 
charge  at  50  ft,  with  Cp  =  2.25.  The  bubble  appears  to  come  to  a 

rather  abrupt  halt  shortly  after  the  minimum.  Although  this  sort 
of  motion  was  not  described  in  early  photographic  studies  of 
bubble  motion,  the  actual  motion  of  the  bubble  'centre'  in  such 
studies  is  very  difficult  to  measure  because  of  the  drastic  shape 
changes  which  are  occurring.  Such  behaviour  has  however  been 
reported  in  some  recent  studies  [07], 

To  further  confirm  the  usefulness  of  the  modified  migration 
equations,  they  were  integrated  for  several  other  charge  depths. 

The  results  are  summarized  in  the  following  table. 
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Charge 

Depth 

(ft) 

Migration 

(ft) 

- *1 

Migration 

P1 

Max  Bub.  Radius 

P 

(0)  i 

1 

l 

Calc  j  Meas 

Calc 

Meas 

Calc 

Meas 

1  50 

26.9  i  24.2 

1.15 

1 .06 

0.18 

0.15 

i  100 

15.7!  15.1 

0.80 

0.78 

0.40 

0.33 

!  150 

10.9  11.0 

0.62 

0.63 

0.56 

0.63 

;  450 

i 

i 

i  3.2  4.2 

1 

0.26 

0.32 

0.92 

0.97 

Again  the  agreement  is  very  reasonable.  The  underestimation  of  the 
migration  for  the  deep  case  is  to  be  expected  since  little  overall 
bubble  distortion  will  occur  in  this  case  and  the  value  of 

will  accordingly  be  too  high. 

Effect  of  Migration  on  Whipping  Response 

As  the  modified  equations  seem  to  give  results  in  reasonable 
accord  with  measured  results  for  the  migration  and  the  far  field 
pressures,  they  are  adopted  here  to  show  the  sort  of  differences 
in  whipping  response  to  be  expected  for  migrating  bubbles.  Figure  04 
shows  the  surface  accelerations  near  a  5001b  charge  of  TNT  at  ?0  ft 
depth  and  shows  a  considerable  change  in  the  bubble  pulse.  Directly 
over  the  charge  the  peak  acceleration  at  the  pulse  is  slightly  reduced 
L>!t  the  duration  of  the  positive  phase  is  greatly  increased.  To  the 
side  Ute  shape  becomes  mere  similar  to  the  non-migrating  pulse  shape 
(aithough  still  of  increased  duration)  but  the  peak  acceleration  is 
much  reduced.  Apart  from  the  slight  change  in  the  phasing  o.  the  bub¬ 
ble  pulse  relative  to  the  ship  motion,  the  only  significant  change  in 
the  forcing  function  is  in  the  bubble  pulse.  Since  the  positive 
phase  of  this  is  still  much  shorter  than  the  basic  ship  period  it  can 
be  considered  impulsive.  The  bubble  pulse  impulse  and  also  the 
initial  expansion  impulse  for  both  a  1001b  and  the  5001b  charge  have 
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been  calculated  for  a  large  selection  of  depths  and  horizontal  stand¬ 
offs.  The  reMilting  impulses  have  been  normalised  to  an  arbitrary 
impulse  and  the  resulting  relative  impulse  contours  are  shown  in 
Figures  05  and  06.  In  comparing  the  initial  impulses  with  the  bubble 
pulse  impulses  it  should  be  remembered  that  for  a  non-migrating  bub¬ 
ble  with  no  energy  losses  the  latter  impulse  is  exactly  twice  the 
former.  At  shallow  depths  the  bubble  equations  break  down  and  the 
position  of  the  contours  shown  by  the  dashed  line  is  conjecture  only. 

The  figures  show  clearly  that  over  a  considerable  depth  range  the 
bubble  impulse  is  greatly  enhanced  by  migration  effects:  for  example 
at  a  depth  of  45  ft  the  bubble  pulse  gives  over  seven  times  the 
impulse  due  to  the  initial  expansion.  At  50  ft  the  ratio  drops  to 
5  to  1 .  At  points  not  directly  above  the  charge  the  enhancement  is 
much  reduced  and  can  fall  well  below  the  factor  2  expected  for  a  non¬ 
migrating  bubble.  For  the  more  strongly  migrating  5001b  charge  the 
maximum  horizontal  standoffs  at  which  given  impulse  levels  are  reached 
during  the  initial  expansion  and  during  the  bubble  pulse  are  almost 
the  same.  For  the  1001b  charge  migration  effects  are  clearly  of  much 
less  importance. 

The  extent  to  which  the  enhancements  in  the  bubble  pulse  impulse 
are  due  to  its  increased  duration  is  emphasised  in  Figures  07  and  08 
which  show  the  relative  peak  accelerations  during  the  initial  expan¬ 
sion  and  the  bubble  pulse.  Only  for  a  very  small  depth  range,  and 
only  directly  over  the  charge,  are  the  peak  accelerations  increased 
by  migration.  At  all  other  points  there  is  a  reduction  in  peak 
acceleration. 

The  impulses  illustrated  in  Figures  05  and  06  are  ''ot  directly 
related  to  the  total  ipimlse  cn  the  ship.  At  any  point  on  the  ship 
axis  the  impulse  per  unit  length  of  ship  will  be  proportional  to  the 
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impulses  shown,  the  constant  of  proportionality  being  the  sura  of  the 
displacement  and  added  mass  (again  per  unit  length)  at  the  point 
concerned.  The  overall  effect  on  the  ship  will  be  an  integral  of 
these  impulses  along  the  ship.  The  contours  do  however  give  some 
idea  of  the  importance  of  including  migration. 

Bubble  Jets 

As  discussed  earlier,  much  of  the  momentum  and  onei  gy  lost  by 
a  migrating  bubble  during  the  bubble  pulse  will  simply  have  been 
transferred  to  a  column  of  water  which  then  will  travel  on  upward, 
initially  at  a  speed  close  to  the  maximum  bubble  vertical  velocity. 

For  charges  directly  beneath  a  ship  such  jets  will  presumably  transfer 
much  of  their  momentum  to  the  ship.  Unfortunately  the  nature  of  such 
jets,  their  length,  diameter  and  transmission  characteristics,  is 
presently  unknown.  Only  one  study  of  any  kind  seems  to  have  been 
made  (reference  07)  and  that  only  considers  the  pressure  pulse  emit¬ 
ted  by  the  jet  at  its  moment  of  creation.  The  pressure  pulse  analysis 
is  rather  dubious  and  since  it  does  not  consider  the  subsequent 
motion  of  the  jet  through  the  water  the  work  is  unsuitable  for  quanti¬ 
fying  the  jet  effect. 

Since  the  motion  of  small  scale  bubbles  in  underpressure  tanks 
is  known  to  closely  resemble  that  of  full  scale  bubbles  (except  for 
the  periods  when  internal  bubble  pressure  is  important),  the  nature 
and  motion  of  any  jets  could  probably  be  investigated  quite  effective¬ 
ly  in  such  tanks.  The  problem  is  not  pursued  in  this  thes?'*. 
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Description  of  Computer  Program  II 


The  general  features  described  in  Sections  L,  N  and  0,  ip.  the 
3-D  virtual  mass,  the  effect  of  a  hole  in  the  target  and  the  effect 
of  bubble  migration,  have  ;d  1  been  incorporated  into  a  computer 
program.  The  program  is  written  in  FORTRAN  and  has  been  implemented 
on  hot'  an  IBM  7090  .md  a  KDF-9  computer.  As  for  the  computer 
program  I,  rhe  program  and  detailed  specifications  for  its  structure, 
me  data  format  anu  units,  and  the  type  of  output  it  gives  are  all 
held  at  the  Naval  Construction  Research  Establishment.  Again  as  for 
the  program  I,  the  urogram  has  been  written  so  that  it  splits 
naturally  into  two  sections,  the  first  evaluating  the  vertical 
vibration  modes  of  the  ship,  and  the  second  using  these  modes  to¬ 
gether  with  additional  details  on  charge  size  and  location,  to 
determine  the  whipping  response.  The  main  differences  from  the 
first  program  are  summarised  below. 

The  clastic  nature  of  the  ship  is  specified  in  the  same  way  as 

for  program  1  but  the  inertial  properties  are  given  in  a  different 

form.  Instead  of  specifying  the  added  mass  and  the  displaced  water 

mass  at  each  node  the  related  Lewis  type  coefficients  C  and  C  are 

defined,  together  with  the  ship  beam.  Since  rotary  inertia  effects 

appear  to  be  generally  small  this  program  does  not  include  them.  The 

program  uses  the  3-D  flow  approximation  equations  to  determine  the 

fluid  inertia  matrix  as  described  in  section  L.  A  Cholcski  routine 

T 

is  used  to  split  the  total  inertia  matrix  into  the  form  LL  and  an 

HGW  eigenvalue  subroutine  then  determines  the  ship  frequencies  and 

i  -T 

transformed  mode  shapes  from  the  symmetric  matrix  L  'S  L  where  S 
is  the  elastic  stiffness  matrix.  This  part  of  the  program  produced 
the  ship  frequency  results  described  in  section  L. 
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Using  these  modes  the  program  will  then  calculate  the  positional 
mode  coefficients  y^(x),  IL(x)  and  S^(x)  described  in  section  1  for 

use  in  computing  the  eventual  response  by  combining  the  individual 
modal  responses.  The  solution  of  the  modal  equations  is  more  com¬ 
plicated  than  for  program  I  since  the  forcing  function  is  no  longer 
proportional  to  the  bubble  volume  acceleration.  Due  to  the  use  of 
the  migrating  bubble  equations  the  forcing  function  is 

£.T  M2  u(t) 

where  M0  is  the  water  inertia  matrix  defined  by  equations  (L20), 

(L22)  and  (L23),  £^is  tl,e  motle  and  u^.(t)  is  given  by 

equations  (013)  and  (01 4) •  In  consequence  a  matrix  multiplication  is 
involved  at  each  time  step  and  for  each  mode.  In  computing  the 

values  of  u(t)  allowance  is  made  (if  so  required)  for  the  effect  of 

a  hole  in  the  target,  as  described  in  Section  N. 

As  before,  the  modal  equations  are  integrated  by  a  standard 
Runge-Kutta  integration  subroutine  and  the  initial  conditions  are 
that  the  mode  coefficients  and  their  first  time  derivatives  are  all 
zero. 
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The  computer  program  II  has  been  used  to  carry  out  an  extensive 
set  of  calculations  and  typical  results  are  shown  in  figures  Ql  to  Q3. 
These  figures  compare  the  results  of  calculations  carried  out  first 
with  migration  included  and  then  without  it.  Figure  Ql  shows  the  effect 
of  migration  on  the  amplitudes  of  the  heaving  mode  and  the  basic 
whipping  mode  for  a  charge  directly  beneath  the  centre  of  the  ship,  and 
also  for  a  similar  charge  at  the  same  depth  but  off  to  one  side. 

Figure  Q2  shows  the  same  charge  and  geometries  but  compares  the  vertical 
midship  displacements  and  the  corresponding  midship  longitudinal  deck 
strains.  Figure  Q3  shows  the  vertical  displacements  and  strains  for 
the  same  charge  at  a  slightly  greater  depth. 

The  most  obvious  feature  of  the  graphs  is  that,  despite  the  great 
enhancement  of  the  local  force  above  the  charge  which  was  shown  in 
section  0  (Figure  06),  the  inclusion  of  migration  does  not  markedly  change 
the  motion  of  the  ship.  At  the  shallower  depth  the  deck  strains  are 
increased  by  2($  but  to  the  side,  and  at  the  greater  depth,  the  motion 
is  little  affected.  This  is  presumably  due  to  the  considerable  length  of 
ship  which  is  affected  by  the  whipping  forces.  Section  0  showed  that 
the  forces  on  the  ship,  although  greatly  increased  above  the  charge, 
were  reduced  at  points  some  way  off  the  vertical.  The  integrated 
effect  of  parts  of  the  ship  away  from  the  charge  evidently  offsets  the 
local  intensified  loading. 

For  deeper  charges,  later  bubble  pulses  will  begin  to  have  an 
effect  on  the  ship  and  in  such  cases  gravity  migration  would  be  much 
more  significant  as  it  could  bring  the  bubble  from  too  far  awt^y  to  be 
important,  in  to  potentially  damaging  ranges.  Such  later  pulses  are 
less  well  known  than  the  first  pulse  and  have  been  excluded  fro*  the 
analysis.  The  bubble  equations  will  predict  several  oscillations  for 
deep  charges  but  the  accuracy  is  not  known. 


Comparison  of  figures  '»l  and  02  shows  that  the  basic  whipping  mode 
is  the  one  principally  responsible  for  the  observed  strains  but  that 
the  higher  modes  still  have  a  significant  effect.  The  higher  mode 
concerned  in  these  calculations  is  the  third  whipping  mode,  since  the 
longitudinal  charge  position  is  close  enough  to  a  node  of  the  second 
whipping  mode  for  it  t'  be  little  excited.  Only  the  first  three 
whipping  modes,  and  the  two  rigid  body  modes,  are  included  in  the 
calculations. 

The  computer  program  II  was  used  to  produce  both  the  migrating 
curves  and  non-migrating  ones  so  that  minor  differences  in  the  bubble 
equations  and  the  ship  representation  should  not  be  included  in  the 
results.  Program  11  has  a  facility  for  suppressing  ihe  migration 
terms  in  the  bubble  equations.  Unfortunately,  since  the  bubble  motion 
is  derived  from  direct  integration  of  the  bubble  equations,  when 
migration  is  suppressed  the  motion  is  that  due  to  an  undamped  bubble 
with  no  energy  losses.  If  such  losses  were  included  the  bubble  pulse 
would  be  slightly  reduced  and  the  amplitudes  for  the  non-migrating  case 
would  be  correspondingly  reduced.  A  second  minor  deficiency  in  the 
graphs  is  in  the  displacement  curves.  The  program  automatically  cuts 
off  the  bubble  pulsations  just  after  the  first  bubble  pulse.  This  is 
because  of  uncertainties  in  the  bubble  forces  after  the  pulse.  In 
consequence  an  upward  impulse  is  imparted  during  the  pulse  but  the 
subsequent  suction  which,  physically,  should  follow,  has  been 
suppressed.  The  shin  then  has  only  gravity  restraining  it  and,  since 
gravity  has  a  much  smaller  effect  than  the  bubble  forces,  the  ship 
displacement  grows  at  a  steady  rate  over  the  time  period  concerned. 

This  is  most  obvious  in  the  heaving  mode  shown  in  figure  Ql(a). 

Since  the  strains  in  the  ship  are  virtually  independent  of  the  rigid 
body  motion  they  are  little  affected. 
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FIGURE  Q  I 


—  MIGRATING  BUBBLE 

-  NON  -  MIGRATING  BUBBLE 
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500  LB.  TNT  45  FEET  DEEP  30  FEET  OFF  CENTRELINE 


COMPARISON  OF  SHIP  RESPONSES  TO  MIGRATING 
AND  NON -MIGRATING  EXPLOSION  BUBBLES 


FIGURE  0  2  b 
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R.  Plastic  Deformation  of  the  Target 

So  far  the  entire  analysis  has  been  based  on  elastic  deformations 
of  the  target.  If  the  stresses  in  the  main  girder  reach  yield  in 
tension  then  plastic  yielding  will  occur  and  the  elastic  equations 
will  no  longer  apply.  Nevertheless,  so  long  as  the  total  duration  of 
such  plastic  yielding  in  each  vibration  cycle  is  short  compared  to 
the  basic  elastic  vibration  period,  the  motion  of  the  ship  will  not 
differ  greatly  from  the  fully  elastic  motion.  If  however  the  stresses 
in  one  of  the  girder  flanges  reach  either  yield  stress  or  the  flange 
buckling  stress,  in  compression,  then  the  flange  will  buckle  and  the 
stiffness  at  the  section  will  fall  dramatically.  Such  failure  has 
been  observed  in  static  tests  on  naval  vessels  (references  [C3]  and 
[C8])  and  the  buckling  has  been  found  to  be  very  localised.  The 
drastic  reduction  in  stiffness  at  the  buckled  section  allows  the 
girder  on  either  side  to  rotate  to  relieve  the  moments  which  caused 
the  buckle.  This  type  of  failure  will  clearly  modify  greatly  all 
subsequent  target  motion.  To  show  the  type  of  motion  which  results, 
approximate  equations  to  describe  the  motion  are  developed  here.  The 
plastic  equations  described  are  not  intended  to  be  particularly 
accurate,  but  merely  to  be  close  enough  to  reality  to  show  the  type 
of  motion  which  is  likely  to  occur. 

With  the  geometry  of  Figure  Rl,  the  equations  of  motion  of  the 
ith  mass  aie  clearly 
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Li  R  i-J 


♦  Pi 


(Rl) 


ri  Yi  =  “Li  ”  **R  i-l  +  Qi 

where  Pi  and  Qi  are  the  external  force  and  moment.  and  are  of 

course  hydrodynamic  and  hydrostatic  in  origin  and  for  simplicity 
the  strip  theory  forces  are  assumed.  Thus 
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and  the  equations  (Rl)  become 
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When  the  deformations  are  elastic,  the  shear  forces  S  and  bending 
moments  If  transmitted  by  the  beans  on  either  side  of  the  lumped  mass 
are  determined  by  the  equations  of  Section  C  and  the  equations  (R3) 
are  identical  to  those  already  considered.  The  plastic  nature  of 
the  problem  is  entirely  determined  by  these  terms. 

In  the  theoretical  model,  since  the  beams  are  considered  to  be 
weightless  the  shear  force  in  each  is  constant  and  the  bending 
moment  is  a  linear  function  of  the  distance  from  the  end.  The 
maximum  bending  moment  will  consequently  necessarily  occur  at  one  of 
the  ends  and  yielding  will  start  at  such  a  point.  In  a  physical 
model  yielding  would  then  spread  outward  along  the  beam  but,  to  keep 
the  model  as  simple  as  possible,  it  is  assumed  here  that  the  dis¬ 
tributed  ’excess*  curvature  near  the  mass  is  equivalent  to  a 
discrete  rotation  of  the  beam  end  relative  to  the  lumped  mass.  The 
rotation  is  assumed  to  take  place  at  a  constant,  plastic,  bending 
moment.  When  the  rotation  exceeds  a  specified  amount  it  is  assumed 
that  the  beam  will  buckle.  This  again  is  assumed  to  be  simply  a 
rotation  of  the  beam  end  relative  to  the  lumped  mass  but  at  a  much 
lower  value  of  the  constant  plastic  moment.  The  complete  moment 
plastic  rotation  diagram  assumed  is  illustrated  in  Figure  R2.  The 
fully  elastic  motions  already  considered  involve  no  plastic  rotations 
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and  so  lie  entirely  on  the  line  GOA.  For  the  plastic  case,  if  the 
lumped  mass  displacements  and  rotations  are  such  that  the  moment  at 
one  end  of  a  beam  section  is  tending  to  increase,  the  moment  will 
rise  elastically  from  zero  at  0  until  it  reaches  the  plastic 
moment  Mpj .  After  this  the  plastic  yielding  will  be  equivalent  to 
a  rotation  of  the  beam  end,  at  constant  plastic  moment  until  the 

point  B  is  reached,  at  a  plastic  rotation  "0^.  Actually  a  more  likely 

route  to  B  is  indicated  by  the  dotted  line,  but  OAB  should  be  a 
reasonable  approximation.  At  B  it  is  assumed  that  the  section 
buckles  and  the  beam  end  instantaneously  rotates  to  reduce  the 
moment  at  the  end  to  The  moment/rotation  curve  ‘jumps*  from 

B  to  B'.  The  rotation  then  increases  at  constant  moment  M  „  until 

pZ 

the  tendency  for  the  moment  to  try  to  increase  further  stops,  spy 
ax,  C.  If  the  tendency  then  is  for  the  moment  to  reduce,  the  beams 
either  side  of  the  hinge  will  unload  elastically  and  normal  elastic 
vibration  could  occur  up  and  down  CD,  A  continued  tendency  to  unbend 
would  lead  to  the  plastic  moment  and  plastic  rotation  along  DE. 

From  E',  the  buckled  flange  would  start  to  straighten  out  and  carry 

an  increasing  load,  possibly  as  indicated  by  the  dotted  line,  but  DEF 

is  a  simple  approximation  for  the  return  to  full  strength.  Plastic 

rotation  would  then  continue  at  the  constant  moment  >11  ,  until 

P« 

the  opposite  flange  buckled,  say  at  H,  and  a  jump  to  the  point  H* 

occurred.  Along  any  of  the  horizontal  branches  of  the  curve,  a 

reversal  in  the  tendency  of  9  to  change  would  lead  to  elastic  motion 

P 

along  a  vertical. 
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The  whole  moment/rotation  relation  is  clearly  highly  over¬ 
simplified;  for  example  the  negative  plastic  aoaent  would  really  be 


of  different  Magnitude  to  the  positive  ones,  the  section  would 
buckle  at  a  lower  moment  the  second  time  than  the  first  since  the 
plastic  strain  reversal  between  the  two  would  not  completely 
straighten  the  flange,  etc.  As  a  simple  model  however  it  probably 
contains  most  of  the  essential  features. 

For  shear,  since  the  shear  force  is  constant  along  each  beam, 
it  is  reasonable  to  assume  that  the  beam  acquires  a  uniform  shear 
angle  along  its  whole  length  once  the  shear  force  exceeds  the  full 
plastic  shear  force.  A  shear  force/shear  angle  relation  similar  to 
that  for  the  moment  has  been  assumed  although,  in  this  case,  the 
reduction  in  stiffness  is  probably  less  dramatic. 

To  produce  an  algorithm  quantifying  the  above  characteristics, 
the  plastic  rotations  ©Li,  0Ri  and  0gi,  as  defined  in  Figure  R3,  are 

introduced.  In  terms  of  these  angles  it  is  possible  to  define  a 
deflection  and  rotation  of  the  right  hand  end  of  the  beam,  relative 
to  the  left,  in  the  form 


»R  *  *i,1  <Yi  *  9Ri  *  9sl)1 

VR  *  (yU1  *  9Li*l'  "  *Yl  *  9Ri' 


These  are  the  deflections  which  give  rise  to  elastic  stresses  within 
the  beams.  The  shear  force  and  bending  moments  in  a  beam  with  these 
relative  deflections  arct¬ 
ic  =  ot  [3yR  +  (c  -  0  «vR] 

UR  =  -  l*L  =  U  C-3yR  +  (e  +  2)  *vR]  (R5) 

SR  =  \  »  6a[yR  -  \  «vR] 


where,  as  in  Section  C, 
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and  SR  are  straight 

lines  of  slopes  -3/(c  -  0#  3/(«  +  2)  and  2  respectively,  so  that 
the  elastic  region  in  the  (yR,  4vR)  plane  has  boundaries  of  the  fora 

shown  in  Figure  R4«  Initially,  yR  and  vR  will  be  zero  and  the 

bending  moments  and  shear  force  will  all  be  zero.  When  the  explosion 
occurs  the  y^  and  will  start  to  change  and,  for  each  beaa  section 

the  values  of  yR  and  4vR  calculated  froa  (R4)  will  describe  a  curve 

in  the  (yRt  «vR)  plane.  The  motion  of  the  ship  will  be  entirely 

elastic  until  one  of  these  curves  meets  a  boundary.  After  this  the 
value  of  the  appropriate  plastic  rotation  (0Ri»  6,j)  will 

begin  to  vary  in  such  a  way  that  the  appropriate  force  or  aoaent 
remains  constant.  The  point  in  the  (yR,  £vR)  plane  will  either 

remain  fixed  at  the  boundary  or  will  slide  along  the  boundary.  In 
any  practical  computing  scheme  the  integration  of  the  equations  of 
motion  will  proceed  in  a  series  of  time  steps  and  at  each  step  the 
corresponding  point  in  the  (yR,  6vR)  plane  will  jump  to  a  new 

neighbouring  position.  Consequently  instead  of  sliding  along  the 
boundary  the  point  nay  jump  just  outside  it.  It  may  then  be  brought 
back  to  the  boundary  by  adjustment  of  the  appropriate  plastic 
rotation.  Changes  in  the  different  rotations  move  the  point  in  the 
directions: 


In  the  (yR,  £vR)  plane  curves  of  constant 
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0Ri  »  at  ^5 °  to  the  yR  axis 
\i+t  »  Parallet  to  the  4vR  axis 

0^  ;  parallel  to  the  yR  axis 

In  any  of  the  regions  narked  I,  only  one  of  the  plastic  conditions 
is  violated  and  there  is  one  way  to  bring  the  point  back  to  the 
elastic  boundary  so  that  its  position  is  uniquely  defined.  In  any 
of  the  regions  narked  II,  two  conditions  are  violated  and  the  eventual 
posit  ion  of  the  point  on  the  boundary  depends  on  the  order  in  which 
the  two  adjustments  are  aade.  The  non-uniqueness  is  not  noraaliy  a 
problem  however  since  the  point  will  always  end  up  close  to  the 
corner,  which  is  where  it  would  have  been  had  its  notion  been 
continuous  rather  than  in  discrete  steps.  The  only  tine  when  non¬ 
uniqueness  is  really  a  difficulty  is  when  buckling  occurs  and  one 
pair  of  boundary  lines  suddenly  moves  in  to  the  new  plastic  noaent 
position  and  leaves  the  point  far  outside.  In  this  case  it  is  clear 
which  rotation  should  be  anended  first  and  if  the  novenent  of  the 
boundaries  is  considered  as  a  continuous  process,  the  eventual 
location  of  the  point  is  obvious.  If  the  anended  position  of  the 
point  takes  it  outside  one  of  the  other  boundaries  it  should  be  left 
at  the  intersection  of  the  two  boundaries  which  is  where  it  would 
have  ended  had  the  transition  been  continuous.  This  particular 
problem  is  principally  a  computer  programming  one.  In  the  regions 
marked  III,  all  the  plastic  conditions  are  violated.  In  general 
such  regions  cannot  be  reached  with  a  reasonable  time  step  length, 
except  for  very  exceptional  boundaries. 

Unlike  solid  beams,  a  ship's  hull  girder  is  usually  designed, 
in  order  to  give  an  efficient  structure,  to  have  nost  of  the  cross- 


sectional  material  as  far  as  possible  from  the  neutral  axis*  In  a 
small  ship  (eg  a  destroyer)  most  of  the  strength  (cross-sectional 
inertia)  is  provided  by  the  main  deck  and  the  bottom  near  the  keel* 
with  the  neutral  axis  about  halfway  between*  In  consequence*  instead 
of  having  a  moment/ curvature  relation  which  changes  slowly  from  fully 
elastic  to  fully  plastic*  as  for  a  solid  beam*  the  ship  behaves  fully 
plastically  almost  as  soon  as  it  starts  yielding  at  the  outermost 
fibres.  It  has  been  found  convenient  in  the  computer  program  to 
define  the  plastic  moments  and  *P2  in  the  form 


M  =  f  M 

pi  i  y 


(R6) 


**P2  3  f2  My 


with 

Io 

H  s  — *  ,  the  moment  at  which  yielding  first  starts*  j/y 

y  y 

being  the  section  modulus*  Similarly*  for  shear  failure 


sPi  =  f3  sy 


Sp2  *  f4  Sy 


(R7) 


withSy  «Agery 


The  factors  f.  to  t .  can  be  specified  as  part  of  the  data  for  the 
1  4 


program. 

One  further  detail  regarding  the  plastic  mcments  is  that*  since 
the  initial  yielding  at  the  moment  really  represents  plastic 

curvature  in  the  region  around  lh«  lumped  mass  instead  of  a  rotation 
at  it  as  assumed  by  the  program*  the  changeover  to  in  the 

program  has  been  arranged  to  occur  then  the  total  rotation 


at  a  aass  exceeds  a  given  critical  rotation^.  In  this  way  the 

program  ensures  that  both  sides  of  a  lumped  mass  reduce  in  plastic 
moment  together. 

Because  the  rotations  of  the  masses  in  the  plastic  case  form  a 
vital  part  of  the  analysis  and  can  no  longer  be  readily  related  to 
displacements  elsewhere,  as  in  the  fully  elastic  case,  it  is  not 
possible  to  halve  the  number  of  equations  by  dropping  all  rotary 
terms.  In  consequence,  instead  of  having  say  five  second  order 
equations  to  solve,  as  in  the  normal  mode  elastic  analysis,  each 
plastic  integration  involves  two  second  order  equations  for  each 
lumped  mass,  ie  forty  equations  for  a  20  mass  model.  This  means  that 
the  time  step  oust  be  very  small,  being  restricted  to  a  small  fraction 
of  the  shortest  period  present  (rotary  oscillations  of  a  single 
lumped  mass  with  a  consequent  very  short  period).  The  combination 
of  a  large  number  of  equations  coupled  with  a  very  short  time  step 
leads  to  excessive  computing  times  and  a  simple  damping  technique  was 
introduced  to  control  the  instability  of  the  integration  routine  and 
enable  the  normal  step  length  stability  criterion  to  be  violated 
safely.  Linear  velocity  damping  was  introduced  between  each  mass  and 
its  two  adjacent  masses.  Such  damping  affects  principally  the 
independent  motions  of  each  mass  while  leaving  the  overall  low 
frequency  modes  largely  unaltered.  The  final  equations  solved  were 

(*i  *  v*  y't  *  s*i  -  su  *  (*.i  *  v>  “i  -  Vi  * 

(,i  *  *.i>  ad  <n-i  *  -  ^i> 

*ri  *  r»i>  Yi  *  "si  "  *4.1  *  <r„i  *  r«i*  “I  *  kriYi  * 

(*8) 


(ri  *  r.i>  <Vd  <n-i  *  yui  -  2n> 


With  this  technique  it  cannot  be  expected  that  motion  equivalent 
to  vibrations  of  higher  modes  will  be  reproduced  but  the  features  of 
interest  in  whipping  should  be  preserved.  The  damping  constants  were 
chosen  so  that  the  predicted  strains  for  elastic  cases  were  similar  to 
those  predicted  by  the  computer  program  II.  Even  with  the  damping 
included  the  largest  incremental  time  consistent  with  stability  was 
found  to  be  considerably  smaller  than  for  the  corresponding  normal  mode 
case  and  the  computing  times  were  found  to  be  about  20  times  as  great. 

Since  the  method  of  solution  of  the  basic  equations  for  the 
plastic  hinging  program  is  very  different  from  that  for  the  elastic 
programs  1  and  II,  it  is  of  interest  to  compare  the  results  given  by 
the  two  types  of  analysis  for  a  case  where  plastic  deformation  does  not 
occur.  The  results  of  such  calculations  are  shown  in  figure  R5.  The 
agreement  is  very  good  for  the  early  motion  but  falls  off  considerably 
later  in  the  motion  due  to  the  damping  included  in  the  plastic  program. 
The  damping  is  actually  a  little  too  severe  but  the  calculation  shows 
the  general  agreement  expected.  It  also  shows  that,  as  expected,  the 
high  modes  are  more  heavily  damped  than  the  low  ones.  As  for  the 
earlier  programs  the  bubble  effects  have  been  suppressed  shortly  after 
the  first  pulse. 

Figure  R6  shows  a  similar  pair  of  calculations  for  a  larger  charge, 
where  the  elastic  limit  has  been  exceeded.  The  strains  are  those  at  the 
middle  of  the  ship,  very  close  to  where  buckling  occurs  in  the  plastic 
calculation  so  that,  after  buckling,  very  little  bending  moment  is 
present  and  the  longitudinal  strains  are  small.  This  would  not  be  true 
at  other  parts  of  the  ship.  The  displacement  curve  shows  that  the 
motion  of  the  centre  of  the  si.ip  is  altered  drastically  by  the  plastic 
hinging.  This  is  even  clearer  in  figure  R7  where  the  ship  profiles  for 
the  same  example  are  shown.  It  is  of  interest  that  in  this  case  the 
charge  geometry  is  such  that  the  hinge  occurs  during  the  suction  phase 
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of  the  bubble  and  bends  towards  the  explosion  rather  than  aw^y  as  might 
be  expected.  The  later  impulse  from  the  bubble  pulse  is  also  very 
evident*  It  occurs  vej"y  close  to  the  ship  and  the  loading  is  fairly 
local.  The  high  shear  stresses  indicated  at  this  time  lead  to  heavy 
local  plastic  shearing  deformations.  Although  the  figures  do  not  show 
it,  the  stresses  were  high  throughout  most  of  the  length  of  the  ship 
and  minor  yiel  ding,  both  in  shear  and  bending,  occurred  at  numerous 
points.  The  plastic  rotation  needed  to  allow  buckling  was  reached 
first  at  the  centre  and  then  the  deformation  there  partially  relieved 
the  stresses  elsewhere. 
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NOMENCLATURE  AND  fGN  CONVENTION  FOR 
THE  PLASTIC  ANALYSIS  IDEALISATION 


FIGURE  R! 


IDEALISED  MOMENT/ PLASTIC  ROTATION  RELATION 

FIGURE  R2 


DEFINITION  OF  PLASTIC  ANGLES 
FIGURE  R3 
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CHARGE  WEIGHT  60  LB.  TNT  AMIDSHIP  AND  45  FEET  DEEP 


FIGURE  RS 


CHARGE  WEIGHT  500  LB. TNT  AMIDSHIP  AND  45  FEET  DEEP 


COMPARISON  OF  RESULTS  FROM  PLASTIC  AND  ELASTIC  WHIPPING 

PROGRAMS:  PLASTIC  MOTION 


Fifcy.fifi-.ftA 
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O  PLASTIC  HINGE  ANGLE  >  I 
////PLASTIC  SHEAR  ANGLE  > '/2° 

SHIP  PROFILES  DURING  PLASTIC  HINGING 
SOO  LB.  TNT  CHARGE  45  FEET  DEEP  BELOW  SHIP  CENTRELINE 


COMPARATIVE  PROFILES  FOR  PLASTIC  AND  ELASTIC  WHIPPING 


SHiP  PROFILES  DURING  ELASTIC  MOTION 
SOO  LB. TNT  CHARGE  <5  FEET  BELOW  SHIP  CENTRELINE 

COMPARATIVE  PROFILES  FOR  PLASTIC  AND  ELASTIC  WHIPPING 

FIGURE  R7b 
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Appendix 

Approximate  Analysis. for  Bilge  Keel  Effect 
Consider  a  ship  section  as  shown  in  Figure  1.  The  free  surface 
can  be  allowed  for  in  the  usual  way  (for  whipping  frequencies)  by 
considering  the  flow  around  the  underwater  section  and  its  image 
(the  dotted  jurve)  in  the  free  surface.  Suppose  that  the  potential  <f>Q 
is  known  for  the  flow  around  the  section  without  the  bilge  keels.  If 
the  length  d  of  the  bilge  keel  is  very  saall  compared  to  the 
dimensions  of  the  section  then  its  presence  will  not  greatly  alter 
the  main  flow  and  the  new  potential  will  be  <t>  =  ^  wher<*  9^  is 

only  significant  close  to  the  keel  and  is  the  disturbance  potential 
necessary  to  adjust  the  original  flow  to  the  flow  around  the  keel. 

In  the  original  flow,  without  the  keel,  the  fluid  at  the  bilge 
keel  position  nad  the  velocity 


u 

n 


-  U  cos  0 

o 


away  from  the  plating  and  the  velocity 
d<4 

u  =  ~  -7—  =  -  all 
s  ds 

tangential  to  it.  The  factor  a  (positive)  can  be  found  from  the 
known  potential  <PQ, 

If  the  length  of  the  bilge  keel  is  small  compared  ir  the  local 
radius  of  curvature  of  the  section,  the  flow  local  to  the  bilge  keel 
attachment  point  will  be  almost  like  the  flow  past  a  flat  plate  as 
in  Figure  2.  The  motion  of  both  water  and  plate  normal  to  the  plate 
is  ignored.  It  is  clear  that  the  relative  velocity  of  the  plate  and 


water  is 


Tne  necessary  fora  of  the  disturbance  potential  <f>  is  then  that  for 

tr.e  flow  along  a  wall  wit!)  a  projecting  plate  as  in  Figure  3»  This 
flow  is  identical  to  the  flow  normal  to  a  flat  plate  of  twice  the 
width  and  is  a  well  known  flow.  The  added  mass  of  a  plate  of  width 
2 

2d  is  xpd  per  unit  length,  so  on  the  half  pertinent  to  our  problea 
the  normal  force  on  the  bilge  keel  will  be  approximately 

^xpd2  (a  +  sin  8^)  U 

This  force  is  perpendicular  to  the  keel  so  that,  resolving  vertically, 
the  downward  force  on  the  section  will  be 

2 

Jfopd  (a  +•  sin  8  )  sin  8^  U  . 

0  0 

In  the  flow  normal  to  a  plate  there  is  also  a  pressure  distribution 
induced  along  the  centre  line  of  the  flow.  In  the  problem  of  flow  past  a 

piato  on  a  wall, along  the  plate  this  pressure  distribution  will  act  on  the  wall 

(ship  bottom)  and  a  component  of  this  pressure  distribution  will  be 

vertical.  However,  neglecting  Bernoulli  pressures,  as  is  customary 

in  added  mass  calculations,  the  positive  pressures  in  front  of  the 

plate  will  be  compensated  by  negative  pressures  behind  and  there 

will  be  no  nett  force,  only  a  couple,  on  the  ship  bottom.  Even  the 

couple  is  balanced  overall  by  that  from  th*  other  bilge  keel.  The 

plate  normal  force  component  can  therefore  be  ignored. 

The  total  downward  force  on  the  whole  section  (image  included) 
due  to  the  bilge  keels  will  be  four  times  that  given  above,  so,  if 
C  is  the  added  mass  coefficient  for  the  section  without  keels  and  C' 
the  coefficient  with  the  keels,  it  follows  from  the  force  equation 
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tr.v'.t 


+  .iicpd  (a  + 


sinO  )sin0 

o 


/V"  (:.Y  2 

~+p\-j  C'  =  r*P^2j  ^  +  (a.  ♦  sinG^JsinO. 

anti  t.ic  increase  in  the  added  mass  coofTicient,  6C  is 

6C  =  J'  -  C  =  2^~J  (a  ♦  sinOo)sin0o 

Unfortunately,  although  this  analysis  is  quite  appealing  physic¬ 
ally,  ii  ..ives  an  order  of  .magnitude  rather  than  a  close  approximation 
to  t.;e  -rue  value  and  some  judgement  is  necessary  to  apply  it.  This 
xs  Joc.rise  for  two  dimensional  flows  the  pressure  near  a  source  (to 
w..ich  the  hilge  keel  is  equivalent)  decays  so  slowly  that  the  integral 
over  the  ncijh  oaring  curved  wall  mqy  have  some  effect.  This  can  be 
snown  *y  considering  the  flow  about  a  circular  cylinder  of  diameter  3 
with  s.'.n.il  fins  of  length  d  at  either  end  of  the  diameter  transverse 
to  t.  e  ;  jtion.  The  .cthod  used  in  section  D  for  rectangles  would  give 
a  modified  added  mass  of  approximately 


U  J 

s  +  d 


so  t.ust  t..c  increase  ov^r  the  usual  value  is,  to  the  lowest  order  in 


\/\ 


*  *! 

:  B 

6y  the  ethod  roposed  were,  when  the  cylinder  moves  upward  with  velo¬ 
city  u  the  water  moves  downward  with  velocity  U  at  the  positions  of  the 
fins,  .’nc  relative  velocity  of  the  fin  and  water  is  therefore  2U  and 
the  total  additional  mass  is  double  t:iit  of  a  [ilate  of  width  2d,  ie 

2*pd2 

*; 

(T 

so  t.wit  ~  «  S-r’,  which  is  of  lower  order  in  d/&  than  before. 


The  exact  solution  to  this  problem  has  been  given  by  Taylor  (D2) 
and  ine  result  is  actually 


O 


d;uble  the  revious  result.  Direct  examination  of  Taylor's  solution 

shows  taut  the  force  on  the  plntc  itself  is  in  fact  equivalent  to 

o  ■> 

Sc“  3",  the  additional  effect  being  due  to  increased  pressure  around 
the  curvature  of  the  cylinder. 

Ti.is  pro’;. era  shows  that  the  order  of  magnitude  given  by  the  above 
method  is  correct  and  that  the  force  on  the  keel  itself  is  likely  to 
he  very  close  to  the  value  calculated  by  the  approximate  method,  but 
t..at  tiicre  :.ay  oc  some  additional  effect  of  a  similar  magnitude, 

Je  ending  on  the  geometry  of  the  hull  near  the  keel. 
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